{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 2c - Monte Carlo Raytracing Methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "nbsphinx-toctree"
    ]
   },
   "source": [
    "### July 2024"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial shows one way in which Optiland can be used to run Monte Carlo simulations. Monte Carlo simulations are a computational technique used to model and analyze complex systems that involve randomness or uncertainty. They rely on repeated random sampling to estimate the behavior of a system and provide insights into the range of possible outcomes.\n",
    "\n",
    "In a Monte Carlo simulation, a large number of random samples are generated based on the input parameters and their associated probability distributions. These samples are then used to simulate the system multiple times, allowing us to observe the distribution of outcomes and calculate various performance metrics.\n",
    "\n",
    "Here, we will investigate how manufacturing variations of a simple singlet lens can impact performance. Note that this is a common tolerancing task.\n",
    "\n",
    "Manufacturing variations to consider:\n",
    "\n",
    "- Lens radius of curvature\n",
    "- Lens index of refraction\n",
    "\n",
    "Performance metrics to consider:\n",
    "\n",
    "- RMS spot size\n",
    "- RMS wavefront error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from optiland import analysis, materials, optic, wavefront"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We choose arbitrary ranges for the lens refractive index and radius:\n",
    "\n",
    "- Refractive index: normal distribution, mean = 1.5, std=1e-3\n",
    "- Radius of curvature: normal distribution, mean = 100, std=0.5\n",
    "\n",
    "In total, we will simulate 1000 systems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_systems = 1000\n",
    "\n",
    "# set random seed for repeatability\n",
    "np.random.seed(42)\n",
    "\n",
    "refractive_index = np.random.normal(loc=1.5, scale=1e-3, size=num_systems)\n",
    "radius_of_curvature = np.random.normal(loc=100, scale=0.5, size=num_systems)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's define a helper class to easily generate random lenses. This is based on the simple singlet in optiland/samples/simple.py."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SingletConfigurable(optic.Optic):\n",
    "    \"\"\"A configurable plano-convex singlet\"\"\"\n",
    "\n",
    "    def __init__(self, radius_of_curvature, refractive_index):\n",
    "        super().__init__()\n",
    "\n",
    "        # define the material for the singlet\n",
    "        ideal_material = materials.IdealMaterial(n=refractive_index, k=0)\n",
    "\n",
    "        # add surfaces\n",
    "        self.add_surface(index=0, radius=np.inf, thickness=np.inf)\n",
    "        self.add_surface(\n",
    "            index=1,\n",
    "            thickness=5,\n",
    "            radius=radius_of_curvature,\n",
    "            is_stop=True,\n",
    "            material=ideal_material,\n",
    "        )\n",
    "        self.add_surface(index=2, thickness=196.667)\n",
    "        self.add_surface(index=3)\n",
    "\n",
    "        # add aperture\n",
    "        self.set_aperture(aperture_type=\"EPD\", value=25.4)\n",
    "\n",
    "        # add field\n",
    "        self.set_field_type(field_type=\"angle\")\n",
    "        self.add_field(y=0)\n",
    "\n",
    "        # add wavelength\n",
    "        self.add_wavelength(value=0.55, is_primary=True)\n",
    "\n",
    "        self.update_paraxial()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's view the first random system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1UAAACVCAYAAACw/9DvAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAANNdJREFUeJztnQmUHGXV92/v0z17lplJyErAxEAIiMAX8/odkHyEyItBXFBcwnJEtuPBgAseZBMJi0cR5QAfCokeP8SFoLLkiJCwaNgCQRREwxvCkmQmM5mlp6f3ru/cp7qqn6qu6ul9m//vnEr1U1U9U5l+urv+de/9X4eiKAoBAAAAAAAAACgKZ3FPAwAAAAAAAADAQFQBAAAAAAAAQAlAVAEAAAAAAABACUBUAQAAAAAAAEAJQFQBAAAAAAAAQAlAVAEAAAAAAABACUBUAQAAAAAAAEAJuEt5cjOSSqVo79691N7eTg6Ho9anAwAAAAAAAKgR3NI3GAzS7Nmzyem0j0dBVJlgQTV37txanwYAAAAAAACgTnj33Xdpzpw5tvshqkxwhIp56392649BY0YcR0ZGqKurK+ddBQBKBXMNVAvMNVAtMNdAtUg1wFzjKNWiQxdOqgsgqkxoKX/8h+vo6Kj16YAS3qTJZFK8hvX6JgXNAeYaqBaYa6BaYK6BapFqoLk2WVkQRFUdc+n9OykcS9b6NBoShRSKxxM0fPAgRaNh8jhS5KYUuR2K+tiRIo8Ypx87FLFfPS5JRy45nFaf9L+pvcVDXnd9v8kBAAAAAEBtgaiqY9p8bnI5YZZRFIpCMadCQ7EQUTJFCU8LhRUHJVIOiitOiou1+jhF2X/j37+o0HUvPiUet3ic1NHioY4WtxBZnX51rY7d1On3iLU4xu82HMvb8RoCAAAAADQ3EFV1zE1nHlnrU2jocPLw8DDd+/MnqXduL/2vFUfbHhtPKhTlJZGiSEKh517cQUNjE7Rm7adoLBynsUiCgpEEjYbjYj0WiVP/WIT+3a+NEzQeTeQUx7kEmDp2U4dfFm7q9lavCy6UAAAAAAB1DkQVaFoURc3T9fv9OY/zuBxiafOqaX693jj53CE6eUlP3r8rmVKEwApGVBGmiTEWYEZBpu57e2gic2wkTpF4yvLncpSrnUWZFB3LiC/r6JgcSWvxuAr8qwEAAAAAgEKBqAJNS2g8JNatra0V/10sfroCHrEUQyyRyhJkQoyxKAurwksXaeEE7RsJ6+OxcIISKcXy53I9mEGApUWZiJppUTJ5LARaJoLmcaGeDAAAAABgMiCqQNNycPhg1URVqbD4md7mE0sxTek40qULr3RUTAgyMc6IMhZuQ+Mx2j0Y0se8Vqw1GQW8LqMAY4Hm9xjGctqiHDXjtEcn6skAAAAAMAVoKFH19NNP06233ko7duygffv20ebNm+mMM84wXFxec801dM899wjP+5UrV9Kdd95Jhx9+eE3PG9SG0ZERsW5ra6Nmhmuu/F6XWHqL6AKQSikUiiUsa8fMgmw0nKC9I2H6VySoR9QmbBwquRRMqyfTIl9yTZlxnBFn2j4WdKgnAwAAAEAj0FCiKhQK0fLly+m8886jM888M2v/LbfcQrfffjtt2rSJFi5cSN/97ndp9erV9Prrr1NLS0tNzhnUtllbo0SqaglHk1SnwuJSFxPJFAWjmYiYWZCZa8oGB0NCnGlRMjYIscItziu3AMuuM1Nryvg5sMIHAAAAQLVoKFG1Zs0asVjBUarbbruNrrrqKlq7dq3Y9otf/IJ6e3vpoYceos997nOWz4tGo2LRGBsb093jeAGNCb92ExMT5Ha7yeF0ivlRzM8Ak8MZfp0iRZA/Tgq/eRGNJ/UomWzuIUfHNBMQFmfvDXM9WSa1MWnz2vrcTkPqolmYyXVmsumHONbnJnee9WQ8T3h+Yb6ASoO5BqoF5hqoFqkGmGv5nltDiapc7N69m/bv30+rVq3St3V2dtIJJ5xA27dvtxVVGzZsoOuuuy5rO6cPsnMcaGznP54DkUikoOd6vV4KBALCkh1UB/4g6nYTdXOmZhs7FvIyeX0ZfxCH4ykajyYpKJaE9DgprO7VtTreOxxV90XUbeM5mmsHvE5q97pECmObj2vLXNTuU8e8Vh+7hGukKxWjnq6QGvHzceqik5xIXQQV+FwLBsdFe3OHA5FYUDkw10C1UBpgrmmZT1NGVLGgYjgyJcNjbZ8VV155Ja1fv94QqZo7dy51dXVRR0cRBSqgbu4qcLooC6tCUz9jsZiIcnV3d1fs/EB9wFb4IU5dNNjfq1ExOXKmRdIGQgn6z2AoHTlL0EQ8aRu9Y3GVMe9IR8A4PTErMiY1kk5H1bjhNOrJgPXdUof4fnI66/PiAzQHmGugWqQaYK65XK6pJaqKxefzicUMv7D1+uKC/EgkEkJQFXtxite/+eGXuMvtoq7Wwl0XNSv8sXCM3usfIocvIKJgVoYf2nggGDLUm3HjaSu4b5qhMbQQYOaxdb0ZizTUkzUv/HmG7ydQDTDXQLVw1Plcy/e8mkZU9fX1iXV/fz/NmjVL387jo48+uoZnBmp598NKMANQLli8TGv1kqO7hbq7Owv6QuDURTbpMPYm02rIMn3JxtLj4YkY7Tk4oVvm8z6b9mTk9zgtTTwmax7N+zi9kfuuAQAAACB/mkZUsdsfC6snnnhCF1Gcyvf888/TRRddVOvTAzUqfITrI6jnO3MtHpdYetqL608WEpExG0EmNZHmY/aPRujf/Zm0Rq43s0O3wpcEWC47fNV1UR23wgofAADAFKShRNX4+Djt2rXLYE6xc+dOmjZtGs2bN48uu+wyuuGGG0RfKs1Sffbs2YZeVmBqoBUVIlIFmhUWLm3cZLnFTbOLtMIXBh4sxlh8pUVYlh0+r8MJ2j04kTk2EhcNp63gKJcstrLEGYsyOTpmqi1jkQkAAAA0Gg0lql566SU66aST9LFmMLFu3TrauHEjffOb3xTmBBdccIFw7/uv//ov2rJlC6IVU5CRtHOfD689AJawZXxXgBcPzS2ynsxSkIlx5rEQZ+mm0XIPs4RN7iKnVJqNO7R6MnmsN4pOG4FoETRPnlb4AAAAwJQVVSeeeGLOfkN85/b6668XC5jajIyOirUfogqAisDiZ3qbTyyFolrhJ/VURCtTD7lX2eB4lHYPhvS0RrbOt/sqCHhdRgGmizDr5tGdUpSM0x65GTYAAADQ1KIKgHwJBtUmzohSAlB/8A2wgJf7ebmpt4jOFSm2wo8lsmrHZHEm73t/JExvRIL6eMKmP5lDt8LPdljMHhsjZrzP70E9GQAATFXyElVyH6d8ueqqq0StEwC1YCIUEusWv7/WpwIAKDMcTRKNlls8dEiXv6h6MtnEI0uQSTVlfNyBA2qUTEt35NRHK9x6PVluy3s5OiYLM1jhAwBAk4uq2267jVasWEFerzevH/rss8/SpZdeClFVIvGk9Rc3mJyx0ISwm863YRsAYGrVk7EVPi/FEI2r/cgy0TCpWbRsj58evzcczph+RBKi6bQVPq4n00WW3Dw6na7IdvepGPVNi1FnwGswAOF9/P8CAABQ5+l/mzdvpp6enryObW9vL+WcQJoVN28TX8CgGDinqIPu+H/vEt/85TvI7ErmcaruZC4HkSe9jffxMdr+sdGZpCS6aM8Dr4qLHHVxibvImbGTfB6XYazuT2/zZLa1aNvSY6QHAdDY8Ht/Ji9FWuFz+qFs4jGW1Zcss21gLEq7BkKqVT7vy2GF3+pzZaJjkomHiI5ZiDTZAIRr0VBPBgAAFRZV9913H3V2dub9Q++++27q7e0t4bQAc+1/f9DWIQvk5umnnqJEMkmHLV5CHPDjvyPfHU4onPrDa4USSaIkr3m7dEw0qFCELfwjCRpKpCiaSIomrZlFHXMKUDxZ+OtjFmcGMWYh2DLHuzJCzWMcG55rEoBsUa2NWUBC1AFQO/j91+pzi2VWZ0vB/fcGhw6SJ9BO47FktqlHOGN5r0XN3jnIUbIxfRy2scJnPaWnLlpFybLGxt5l/JmEzxYAwFQmL1HFluWFcPbZZxd7PkDiv4+aVetTaFj6nxsRqX8fPayt4C/6Z5/5l7DkX7/u05MeyyIsJoktfhzRxnF1bBZk8raIxTGxtGgLRROZ7fGMsFOPT4rfk8MM0/bCyS7y5nVzI1rj2E788XFWYtBrEcGTxR9HAwEAxcPvIRYz3a3F9eDjzw9u/CzXjlk1i9b29Y9FDfVmdjeSPC6HhYmHtamHVZ0Z6skAAFPa/Y+b8fKdM5mOjiKsnAAoM/F4PO8awFIvcPxel1iIPFRNOI2Io2sZoZUWW/GM8DKKNSnSFreOvGnbOD1peCKub4uYxKB2fKHwhZdV5E0INW0sCTvLtEs+bpJjZEGoRfX4d+NOOpjq8Htlmru4ejL+zOH3vtnyflSkLJrFWZyGJ2K0Z0gz+VBFmV3yhd/jtGgWbTT1MEfHtGPZCh83bAAADSeqdu/eLUwotm3bRpEIJ0llPmz5giWZtLaqBaCaJBKJpjep4PcbCwXR7LS4m9YlwbbWbKZijMSZRZu1GNPEnzyWj+G75IZt4lhjtK6Y1NicaZdSSqVd2qV2HAtAbex1OSgWnqDpY0QtXndWhE5bYCIAmuEzh9OJeSnWCl+tJ8sWYLIw00w99o9G6N/9aXfGSJxCUfvrCxZWOQWYpR2+Om71wgofAFADUfXFL35RCKh7771X1E3hgwjUI8lUipxOXMRWEi5q9zk5alQb8cq22LG0qNMjallpkhaROEnwWaZmxlPiIk5OxdRSOrVxtIjUS76TblUDp0XVWrRxLgOUfExSpG188asKPydMCEDN4TnYxk2WW9w0u8j3/Hg0I8q4hkw2+JANP9gAZPfgRObYSFzcyLF7bxos73ktOStmTD+Mfcm0Y2v1GQgAaHBR9eqrr9KOHTto8eLFlTkjAMqAAlHV9HDkh5dA5bM8s+AbS1xbwiIsHEvQwNAwtbS2EV+zyfVvdpG4WA4ByLV0w6FcdXfJogxSOKppHYkzp1TmNkDJ1N1ZO2IaTFLSxyD1EpQDfr93BXgpLtWa30eyIJNNPViEaZb3mkjbOxI2jO2i4zzHdbGVFmL2NWSqQBNCLX2syDYAAEw9UXXcccfRu+++C1EF6hqu9Wv29D9QO1ggeN1qfVjA6yRn3Evd3a1VE/KcRqVF6QwmKTnSLsU2qZbOsk4unqSRiZhR/JnSLlnkFZp5yXrKzgDF0HoglwFKenvGJMVOAGa7ZKLeBjA8J2a0+cRSzI2UcDzjuGhl6iH3Khscj9LuwZBef8ZW+HbRbbazlwVYLlMPud6M15z2iCg0AA0qqn72s5/RhRdeSO+//z4deeSR5PEY7xgdddRR5Tw/AIqCvwAhqkCzwhdRLU61tqUWcBqWVZSN0yQzws0+FdOcgplJ31RNECZzyCwUbiVgJ7zMbQdyGqCk6+km602nu2SiN13TwK9hwOsWS29HYVb42o2QUCxhY+qRncL4/kiY3ggH9THXolmfF4nGz1a1Y4bomEmcaY/9HtSTAVAzUXXgwAF666236Nxzz9W38RsSRhWgnuD5iPQ/ACqbelmkq3dZUi/zS6k0mp3I7Qisjqnn3nSaKUoyFqHujhi1eF15NSdHb7r6uRHC0SVeiPwFP59NgYI20THZ8l4bHwhG9foyHtvNW54fVgKMhZeHEtTTNUKdAa8hOiZH0mCFD0AJouq8886jY445hu6//34YVYC6BqIKgOZOvWyvwe+36k2nm5nk0ZvO3HfOrjddzMIhk8eFSrrJetNl2g7k7k1nrrkzHGNhnKLV0yH1sjxw3RXb4Bdjhc/wnDObemiuilxPpu9LR9DeGw7TcChKodhBkbrI894Knju62JJqxrJNPTTDj0zdGUfY4IoKprSo2rNnD/3xj3+kww47rDJnBECJaL3TIKoAAM3Um44zQQaHhinQ3qGaomSlVE7em06uk8vVm87KIbPcveny7TuXj0mKWRBySicMUjKw8J3JS7sv7+/R4eFh6u7uFn/DUEytJ9OcFc2mHnJN2cBYlHYNhAwW+Xa0+riezGMSZLLpRyal0Rwl4+fi9QUNLao+9rGPCQdAiCpQr8RiMbFGTRUAoJngC0i3y0GtwpzAWTe96eSoWq7edFrdXD696azSOcvdm06uqcuVdmnuYWc4xiQIM5G65ulNx/OODTF4mdVZeD0ZR7nGo/amHuZeZe8cZNfFMX0ctrHC5yCo3Ag6u1l07ubRLNIhykBNRdXpp59OX//61+m1116jZcuWZRlVfOITnyjn+QFQtKgCAADQ3L3p5P505exNp6V0VrI3nW5mUmDfuVzH6A6ZddSbjv8OLGZ4oe7Cn8+vEYsyuXZMNvUw15XtG40YDEDs6sk4kmnXENrguGhhAMICE/VkoGRRxc5/zPXXX5+1D0YVoB5IxONi7USkCgAAmoZ66U2XJc7y7E1nJwDNvemsWiRUpDedKQon7xemKPEYdbUHdVOUljx602k/u5ypl/zzp7mLqyfj143/jtmmHsbomDY+GIrRnqFQpol0JGHbQsKfrieTo2Fmy/vM2CjOWJSh3rD5cBdbrwJAvZJMqcIeH1cAAADKbZDSVvilU8V602kplfn0prOsk7PoTaeJv0g8QfEkiRq8SvSmM/ads+9NZ2WUkqtFgk9KveTXTauB7O0o7u/O9YayiYds+KHb46fX+0Yi9Ka0LRS1DzRokbDJBFjG5CMzbvWinqweqf4nAwAVJpFIiyoYVQAAAGgCqt2bTjaq4Po9u950dumVVqmYmvibrDedVfSv3L3prFoZ2PWmE+JPEnscZZre6rXsTactmuDhvxu7J+rRsXA8Y/IhpTBqBiDcMFodq4KN/z5WcJTLkJ7IYkwy8VDHaTFmEmm8r1YpvM1OUaLqxRdfpK1bt9LAwEBW5OqHP/wh1Zo77riDbr31Vtq/fz8tX76cfvKTn9Dxxx9f69MCVUJLQcU9HAAAAKB5e9Ox+Irk0ZsuUyM3eW86s0NmJXrTmU1SdKFmEn+cIci1fClFEUYtbBSjpaGy4OLzDMeSIprGEcf3h8NCxHFEjYWbnbkL/3yjy2KuZtHGqBkvbPEPyiCqbrzxRrrqqqto8eLFWX2q6iEU+cADD9D69evprrvuohNOOIFuu+02Wr16Nb355pvU09NT69MDVQApqgAAAEBzUK+96bQ0ybx600ltDWRDFK7hMm+THTIjRRikaL3pOEWQ/2ZCEDsdxDrI6XCKG87a5TprrpFwnIZCsbQRjNHh0w5O29QcIbXUxU6/h7p4CWQMP8x1Zlx/1upl99La64W6EFU//vGP6d5776VzzjmH6hGOlH3lK1+hc889V4xZXD3yyCPinL/97W9nHR+NRsWiMTY2pl+Y1/rinO8+/M9gqKbn0Ij0HwjTYCpA744TxYeiBYv9/VEPBZN+eu29kYqdI2ge+C5qMBii9gnY84LKgrkGqgXmWn6woYfX5aI2n6tir0OSI3UJtaZOj1TxmoWYJoRERC0jitT90mNeG8aZnyG2CYHlIBf/vhQLMAexuaGdrhJRsniMBseLc1t2pdtDsKGJ28miTyGP2yXEn9ORXpxE9593nEhrrDX56oGCRRXn1q5cuZLq1Up7x44ddOWVVxrOd9WqVbR9+3bL52zYsIGuu+66rO0jIyM1dzJ8c2CCzr3/jZqeQ+OylP70KhG92l/Ec2eJfx/8vy+U/awAAAAAAKYySRZvCU7pFKP01uwm0a/u3kdHzKpFfNJIMBisjKjiHlVcs8RpdfXG4OCgEEKclijD43/961+Wz2EBxumCcqRq7ty51NXVRR0dRVjFlJHlrR30uwtqP5kajf7+ftqy5TFasngxzZo9u+C7bDtfeYWC4+N09tlfqNg5gma7oxuk9vZ23NEFFQVzDVQLzLVqwbVSmR5sHE3iNEA9kpSOQmWiS+o+Q5RKf26KEunn8ZJIR7fUYzglUR0nxHPTj4toqM1wRZU2LfgnTPZj3E5e0n3Z0m6OXE/mF4uTvA6FOttaRIogR/3afR6x/ZhFs6jVV/tIlSvPFj0Fi6orrriCTjvtNFq0aBEtXbo0q/nvgw8+SI2Ez+cTixmOcNWiY71Ma4uTls3pquk5NCKdqTF6yTlB89odNH+6r/D0P1+cWsJh/O1BAS5ZKeru7qz5ZwZobjDXQLWYSnPNqmZKM7+ITlIzZWdnLzeNtnI0lA0zymlXr9VPccmS2+0kj9tJAVJIUdyUTPH/VRVSWhqhqNtK/x/CsZQQSPY9ubQaKS91Sm6DolbK7D7oz1jDT9aTy+w0WY/ke14Fi6qvfe1rwvnvpJNOounTp9fVHYwZM2YINcmRChke9/X11ey8QHWp1zclAAAAACZ391PNGiZpqhw3OvVZHZPp06WaPlhZvZe7sTJHYLRxd8BrbeWetm9nwwd5zLVGLHziaaGnizfJ5W+cLdp5kSzZhydiwu3P7v/C5yu7+HEj5axeWLrbX3ajYj43UAFRtWnTJvr9738volX1htfrpWOPPZaeeOIJOuOMM3QFzONLL7201qcHqhym5Q9qAAAAAOTG3IdqIpagwYNh8k04RQPgXH2oMm512X2o8o3sFAoHPmTxIkdutJ5SwgGvxU3Tclibq72mrBoPu2z7WXldnMnksI16sbgxNAVOix/ezj25hid4X1gcM2pqIGzXl4p/nZWjXl9ni6H/lOawZ24oLPfOAnUkqqZNmyZS/+oVro9at24dffjDHxa9qbj2KxQK6W6AoPlxu9VpDVEFAACgEUhxZCJpb8FtZ9PNkRg9Vc0qDc3URNey11MiJcRAodgJEY6+yGO+qM+IH5dlnybLRr16FCdb5HCaWyXg64ZQLCkJoJgqeMKq8AnaCCWtuS9Hkexo9bmymvEumB6wtB6Xx3wsPxeiqAlF1bXXXkvXXHMN3XfffRQIBKjeOOuss+jAgQN09dVXi+a/Rx99NG3ZsiXLvAI0LxBVAAAAiklBs2wWG588Dc1OrJgjN9zXyCqSU2wKmq1YkYRJZ8BjSDEzp6Gp0RpjGprXSRQNh2hGdxe1eNVIh5yqxr+7Xi/y+W+sCqDsSJAshoxiSd3HaXV2ApP//3r9UFr49HT46PCetqwmuUIc+TPHsvlCpYQgaGBRdfvtt9Nbb70lRMqCBQuyjCpefvllqjWc6od0v6mLJqpq3WcMAABA/vDFrCxCtFSxyYQKixquObE1FDBFfOwMBQq9D8eaosVOrEjpYtyEVautMUd2NKFilXKmuqSZ0tAk8ZSr+L9UVPMAB3V3t1e9TplNFPSIEAuecIJGNSGUjgrJY00waUKJX3Mr2MDBUEPUojaqnTdNjRaJKJLYbxRG2hh1RaDsokqrVQKgXtHcHFOIVAEAQEHRGlmYGISKVPCf7XSWicBYRnKkbVb1N9r+Yuyds2tq0mJFTinzuKgtXVdiFdlRhY1VGppRMJlT1fgivV6jNbVOZeQ0ODk6JKJGQiBl0uiEMDKPIwlhxmAF/6nbfXKanCp2ejpa9HHGbCFjvqDtY/tuvF6grkQVp/4B0BDpf4hUAQAaTNSwsLATKuFYgoZGxsjbn6CYyQZ6UrMAKe0skkMcFQoLC6t6mExKmSpK+KJ2hluL1hgjO9aGA9Y/Sz7Gk8MwAJRYVxRN0EAwRgdi4zQe5XQ6SfzoNUTmmiL1GBZUdvc0A16XyXbbQ3O6/Yb6oSwnuvS+Vq8brzdoLlEFQL2jpSpwI2gAACilZ425z4xdipmaqpZnz5ocNTrl6lnTYhr7va50bU12ipmovzGkodlFbtR0Na/kgoY6kfqE59SoXC9krh8ymS+YhZJd1JBfe7mmiAVPT7uPFs1szbLoNo95YSEMwJQWVez49+9//1v0gcqHefPm0TPPPEPz588v9fwAKAoO8aOmCoDG71mjiZV8etbYpZyZe9ZofW0q3bPGnC5m27PGom7G6hiPkygSGqee6d3k96qOavVsGABKszgXvYhyRobshZKdTTnXYen1Q+lIENcNHdLVogslUVfkc5MjEaFZM7qoK+BLR5HcYh4CAEoQVSMjI/TYY49RZ2dnPofT0NAQogSg5kBUAVB6z5qs2hcpImM2EMh2NrN2OpvMKrpcPWvMURm7njVaXxtvPj1rTMfk6llTEfMAV4y627xoct4AdUWhmCp8zIYKmgNdrqiRXV0Ro4shKVq0cEarpQOdiBJJY06/y0eEq0YVw9Td3YW5BkC50/+49xMAjQJ/CUBUgabpWWPRa8YqDW0ysSL3rFHtoivbs8bK2cyuZ429WYCxZ41uFV2FnjVgakdMOUqqiZ9s1znJfMFCKPExdm8jv8eZ1ZNodpeflmjNWyXXObP5QpvPXVHXPwBAhUUVLk5Bo8F34hAtBeXqWaM117TqWROJJ2h4dJzcvhBFE8Y+N7msoLPd0srfs8ZKmNj1rGGxolo4T56GZv7ZSEED9Qi/r6wiQeZokUirM48jCdv3Ic93s6HCtFYvLZgRkOqJjM1btWNZFMGaG4DmBEYVoClxOl0QVU3Us0a2c84lVKxT0Ox61pijPyX2rEk3hzRaMBvNAnL1rNHqb/LpWZOxi65OzxoAavUZwOLGrnmrXFtk1eiVI01W8FvFSvz0darW3OYokmzhzcKI33O4iQAAMANRBZoSlwvpf+XsWZNlwZzD6SyXs5m8za7+ptw9a7LSxXL0rDHXzNj1rJFT1dg8YCI4RjOmd5PLhSJuAAzW3LGkJIByR43MQomtue1o9bkkswW1bmjBdLWJqyaIOtOW3eaoEdcVwZobAFBuIKpA0/aqavRIldyzJl+hkrdZQB41OmXrWWORTiZ61rQZe9Zo9TdWKWfZPWvMKWi161nD4j02gfQ30Jzw54ldJEgWQ1lOdHxsNGFbn8fvd13wpNc9HT46rKfN0MdI26fVGKl1RS7U0QEAGldU7d27l2bPnl3ZswGgTHg8nrKIqmJ61kwWyZFTzNSGnZXtWaM7m+XRs8bcbNOuZ012qhp61gBQj8TZmltqzMq1QxmzBbP5grGmiB/b3VzhGyjmnkRdAQ/Nm6ZGi0QUSezPNl/gY1FXBACYsqLqiCOOoDvuuIPOPvvsyp4R0Lnu4TfEXUJQOLtGeymeTNILzw4S1xpzxIeXZIoonl7r23h/UqGEom4PRw+hROoQuve6x8vSsybLgjlHzxo7W+h8+9rAMACA5nOC5DQ4FjgjE1HaeyBIyv44BaNJCso23MJswTSOJGytufljot2X7S7X06HWFVk1b+2QjvV78rPmBgCAqULeour73/8+ffWrX6XNmzfT3XffLRoCg8qyZ2iCJmL2OeXAnoNJLyUTCXJHkuKOKhfxs+DwexzkdvBYvdOqLupjbdt77+yhWGSCTln1sZw9a1QTgtr1rAEANEYaLwsbK9e53I1d1WNYUNmZpnBtkJwqxxGhOd1+SzEkmy9wrVGr143PKgAAqIWouvjii2nNmjV0/vnn09KlS+mee+6h008/vZznAkzcu+7YWp9Cw/LQ5s303nvv0qqPHU7OAs0Dnh19jUaUEfrC8fMqdn4AgMaBU+BY7Fg3a7Vu3ioLJTvTFb7Rw5EfvW6I64rafbRoZmuWZbcY+1ykxEI0p2c6dQS8oo4QAABAAxpVLFy4kJ588kn66U9/SmeeeSZ98IMfFIYAMi+//HK5zxGAgvH7/WIdTyTIB0c2AKY0Ca4rik4WGbIXSlzfaAVHwPVIkGSoMLurRRqn64sszBc4jbdQU5Th4SR1t3pFg3MAAAAN7P63Z88eevDBB6m7u5vWrl2bJaoAqAfa29vEOhIOk8/nq/XpAABKrCsKxbR+REZDBWM6nXXUyK6uiJEtuLU0uYUzWi1rijLmC+qYe46hrggAAABTkCLilL/LL7+cVq1aRf/85z9p5syZ+CuCuqS9vUOsI5EIddb6ZACY4nBdETtdauIn23Uubdkt2XHLQomPsXPC9KetueWeRLM6/bSkz7p5qyyU2nxuNE0GAABQXVF16qmn0gsvvCBS/7785S+X57cDUCG6urt0UQUAKE9d0WTNW3U7bvM4krB10uS6IrOhwrRWL82fHpDqiSyiRn5VFMGaGwAAQEOJKu758/e//53mzJlT2TMCoAxo7pTRaLTWpwJAXcA911jc2DVvlWuLrBq9cqTJCg706BEgyV2ur1Oz5jZGkcxRI3bRRAodAACAKSOqHn/88cqeCQBlxO8PiHUYkSrQRCl0obQ1tyqAckeNzEKJrbntaPWxNXemMSvXDc2fztbcHTmFER/Ltt6w5gYAADDVaRiXCe6T9cgjj9DOnTvJ6/XSyMhI1jHvvPMOXXTRRbR161Zqa2ujdevW0YYNG2CmMUXhu98xRKpAHcHNvO0iQbIYynKi42OjCRFtsoKjPbrg0Zu4+uiwnjZTHyPNbCETLeIGsG5YcwMAAAAl0TBqIxaL0Wc+8xlasWIF/fznP7dMTzzttNOor6+P/va3v9G+fftE7ZfH46Ebb7yxJucMaovL5UL6Hygrcbbmlhqzcu3QyESM9h8cpaRzlMajScl8wVhTxI+5LskKbjpt7knEj+dNC1g2bzWMfYVbcwMAAABgioqq6667Tqw3btxouf/Pf/4zvf766/SXv/yFent76eijj6bvfe979K1vfYuuvfZaEd2ygi+65QvvsbExvR8IL6Ax4deOBXUkEhVpU8X+DNB81tycBmdwnpPc6HQb7vQ+OaWOH0/Era25OflNRH2kSBE/XjQjIMRPlm23bragbvOXWFeEuTo14NeZP8/weoNKg7kGqkWqAeZavufWMKJqMrZv307Lli0Tgkpj9erVIh2Q7d+POeYYy+dxeqAm2GQ4vZCjX6AxUZQUtba2itewUAdAFuCBQICGh4crdn6gOPiDNxxPiYhQUCwJ/bG6TuiPzeNgJClqkuwktrDm9rmozecS0R9+PDPgooXdAX2s7nMJQaQ9bvU6KBUNU2dHOzkc+abR8VnEiZJxiobChHgqyPdzLRgcF/Mn/7kGQOFgroFqoTTAXAsGg1NLVO3fv98gqBhtzPvsuPLKK2n9+vWGSNXcuXOpq6uLOjrUXkeg8eC7ClxLNzo6Si0tLQWnmk5MTIgG16D8cAqc2WlOTZOTm7cat2spdLwkUvbW3J167ZAaEerr8tEHdFe6TBqdwaabBVKLmzxF1hXxXOObMPyZ4XTW5xcCaA7Uu6UOzDVQcTDXQLVINcBc43KSuhdV3/72t+nmm2/Oecwbb7xBS5Ysqdg5+Hw+sZjhF7ZeX1yQH4GAGqnipRizErz+1iS4rkhLodPFUbbbnMFoQRpHbeqKuAmrnCbXyWu/hw7p9ksGCxkhlDFfUI+tZV0Rp+7hMwNUA8w1UC0w10C1cNT5XMv3vGoqqi6//HI655xzch5z6KGH5vWz2KCCmxPL9Pf36/vA1KO9vV2sx8fHxR0QkKkrCsUy9UOyoYLcxNXKsptF00TMPi2Wm7HK1tsshhbOaM30MZKat2bMFtRxq9eFfkUAAAAAaEhqKqpmzpwplnLAroBsuz4wMEA9PT16by1O4Vu6dGlZfgdoLLrS6XvjwWBTiSquK+JGrJr40UwWMq5zactuyY5bFkp8jE0GnVpXZHKXm9XppyV9RqEku9RpQokFFUebAAAAAACmGg1TU8U9qA4ePCjWnM7F/aqYww47TPSkOuWUU4R4+tKXvkS33HKLqKO66qqr6JJLLrFM7wPNz/Tp08Q6FApRPdYVTda8Va8lMo8jCYon7euKzJGgaa1emj89kNW01Rw1YtHkdddn6B0AAAAAoJ5pGFF19dVX06ZNm/Sx5ubHjX5PPPFEUUT28MMPC7c/jlqx8xs3/73++utreNaglvh8LSKdjE0nyg03YVWNE6ybt5qNGMyNXjnSZAUHevQIkGbB7XdTX2dLptbIbx814iawSKEDAAAAAKguDSOquD+VXY8qjfnz59Ojjz5atXMC9Q8bVGiiiqM7UV4SKYoklPSSomiCt6mPedsbY+00FvHSjY/9yyCU5KgS9zqyo9XnyvQhSpsrzJ/OZgsdqvmCFCWS3ej42IDXRU6k0AEAAAAANBQNI6qmIv/ntmdyXryDHChEKUWhaOQISo4R3fz229wBIa8nuqmDvI4Ujewa0qNDPR0+OqynzeQ6pz3ORI24n5G7SGtuAAAAAADQmEBU1TGfP26uqL0BRTaJDYfpwIEBCo6NktehkMepkNeZXosxZbY51LXbwdaeRMcdfwItWrSo1v8NAAAAAADQAEBU1THnrVxQ61No6GZyw8PD1N29tG77HgAAAAAAgOYAV5sAAAAAAAAAUAKIVFmkjTHBYLDWpwJKjFTxa8iukIhUgUqCuQaqBeYaqBaYa6BaNMJc0zSBphHsgKiy+cMtOnRhrU8FAAAAAAAAUCcaobOz03a/Q5lMdk1Bxbx3715qb29Hv58GZmxsjObOnUvvvvsudXR01Pp0QBODuQaqBeYaqBaYa6BajDXAXGOpxIJq9uzZOaNpiFSZ4D/WnDlzan0aoEzwG7Re36SgucBcA9UCcw1UC8w1UC066nyu5YpQadRn8iIAAAAAAAAANAgQVQAAAAAAAABQAhBVoCnx+Xx0zTXXiDUAlQRzDVQLzDVQLTDXQLXwNdFcg1EFAAAAAAAAAJQAIlUAAAAAAAAAUAIQVQAAAAAAAABQAhBVAAAAAAAAAFACEFUAAAAAAAAAUAIQVaDpuOOOO2jBggXU0tJCJ5xwAr3wwgu1PiXQ4Fx77bXkcDgMy5IlS/T9kUiELrnkEpo+fTq1tbXRpz71Kerv76/pOYPG4Omnn6bTTz+dZs+eLebVQw89ZNjPXlJXX301zZo1i/x+P61atYr+85//GI45ePAgfeELXxCNM7u6uuj888+n8fHxKv9PQKPPtXPOOSfrc+7UU081HIO5BvJhw4YNdNxxx1F7ezv19PTQGWecQW+++abhmHy+N9955x067bTTKBAIiJ/zjW98gxKJBNUrEFWgqXjggQdo/fr1wp7z5ZdfpuXLl9Pq1atpYGCg1qcGGpwjjjiC9u3bpy/PPvusvu/rX/86/elPf6Lf/va39NRTT9HevXvpzDPPrOn5gsYgFAqJzym+GWTFLbfcQrfffjvddddd9Pzzz1Nra6v4TOMLEg2+yP3nP/9Jjz/+OD388MPi4vmCCy6o4v8CNMNcY1hEyZ9z999/v2E/5hrIh6eeekoIpueee07MlXg8TqeccoqYg/l+byaTSSGoYrEY/e1vf6NNmzbRxo0bxU2muoUt1QFoFo4//njlkksu0cfJZFKZPXu2smHDhpqeF2hsrrnmGmX58uWW+0ZGRhSPx6P89re/1be98cYb3KpC2b59exXPEjQ6PGc2b96sj1OplNLX16fceuuthvnm8/mU+++/X4xff/118bwXX3xRP+axxx5THA6H8v7771f5fwAada4x69atU9auXWv7HMw1UCwDAwNi7jz11FN5f28++uijitPpVPbv368fc+eddyodHR1KNBpV6hFEqkDTwHczduzYIdJjNJxOpxhv3769pucGGh9OueK0mUMPPVTcreW0BIbnHN+Fk+cdpwbOmzcP8w6UxO7du2n//v2GudXZ2SnSmrW5xWtOw/rwhz+sH8PH82cfR7YAKIRt27aJNKvFixfTRRddRENDQ/o+zDVQLKOjo2I9bdq0vL83eb1s2TLq7e3Vj+Eo/djYmIiW1iMQVaBpGBwcFOFi+Q3I8JgvTAAoFr6I5bSDLVu20J133ikudj/60Y9SMBgUc8vr9YqLDRnMO1Aq2vzJ9ZnGa74IlnG73eLiBfMPFAKn/v3iF7+gJ554gm6++WaRkrVmzRrxvcpgroFiSKVSdNlll9HKlSvpyCOPFNvy+d7ktdVnn7avHnHX+gQAAKDe4QsLjaOOOkqIrPnz59NvfvMbYR4AAACNzuc+9zn9MUcI+LNu0aJFInp18skn1/TcQONyySWX0D/+8Q9DHXKzgkgVaBpmzJhBLpcryz2Gx319fTU7L9B88N21D3zgA7Rr1y4xtzj1dGRkxHAM5h0oFW3+5PpM47XZiIfdsdilDfMPlAKnOvP3Kn/OMZhroFAuvfRSYWiydetWmjNnjr49n+9NXlt99mn76hGIKtA0cCj52GOPFakLctiZxytWrKjpuYHmgi2E33rrLWFzzXPO4/EY5h1bx3LNFeYdKIWFCxeKiwd5bnE9AdevaHOL13xhwjUKGk8++aT47OOIKgDF8t5774maKv6cYzDXQL4oiiIE1ebNm8Uc4c8ymXy+N3n92muvGYQ8Owmynf/SpUupLqm1UwYA5eTXv/61cMbauHGjcCq64IILlK6uLoN7DACFcvnllyvbtm1Tdu/erfz1r39VVq1apcyYMUM4GjEXXnihMm/ePOXJJ59UXnrpJWXFihViAWAygsGg8sorr4iFv5J/+MMfisd79uwR+2+66SbxGfaHP/xB+fvf/y7c2RYuXKiEw2H9Z5x66qnKMcccozz//PPKs88+qxx++OHK5z//+Rr+r0CjzTXed8UVVwjnNf6c+8tf/qJ86EMfEnMpEonoPwNzDeTDRRddpHR2dorvzX379unLxMSEfsxk35uJREI58sgjlVNOOUXZuXOnsmXLFmXmzJnKlVdeqdQrEFWg6fjJT34i3qher1dYrD/33HO1PiXQ4Jx11lnKrFmzxJw65JBDxHjXrl36fr7Avfjii5Xu7m4lEAgon/zkJ8UXCACTsXXrVnGBa17Y3lqzVf/ud7+r9Pb2ihtGJ598svLmm28afsbQ0JC4sG1raxN2w+eee664SAYg37nGF7t88coXrWx1PX/+fOUrX/lK1g1JzDWQD2Qxz3i57777CvrefPvtt5U1a9Yofr9f3MjkG5zxeFypVxz8T62jZQAAAAAAAADQqKCmCgAAAAAAAABKAKIKAAAAAAAAAEoAogoAAAAAAAAASgCiCgAAAAAAAABKAKIKAAAAAAAAAEoAogoAAAAAAAAASgCiCgAAAAAAAABKAKIKAAAAAAAAAEoAogoAAMCUYcGCBeRwOMQyMjJS9d+/bds2/fefccYZVf/9AAAAKgNEFQAAgIZCFiZWy0knnZTz+ddffz3t27ePOjs7qdp85CMfEb/7s5/9bNV/NwAAgMrhruDPBgAAAComTMz88Y9/pAsvvJAuvvjinM9vb2+nvr4+qgVer1f8br/fT9FotCbnAAAAoPwgUgUAAKCh0ISJvAwPD9MVV1xB3/nOd+gzn/lMQT9v48aN1NXVRQ8//DAtXryYAoEAffrTn6aJiQnatGmTSBns7u6mr33ta5RMJvXn8fYbbriBvvzlL1NbWxvNnz9fCLsDBw7Q2rVrxbajjjqKXnrppQr8FQAAANQTEFUAAAAaGq6NYhFz4okn0ve+972ifgYLqNtvv51+/etf05YtW0SK4Sc/+Ul69NFHxfLLX/6S7r77bvrd735neN6PfvQjWrlyJb3yyit02mmn0Ze+9CUhsr74xS/Syy+/TIsWLRJjRVHK9L8FAABQjyD9DwAAQMOSSqXo7LPPJrfbTb/61a9ETVUxxONxuvPOO4UIYjhSxUKqv79fRJyWLl0qarW2bt1KZ511lv68j3/84/TVr35VPL766qvFzzjuuOP0aNm3vvUtWrFihfg5tUo5BAAAUHkgqgAAADQsnO63fft2euGFF0StVLFwyp8mqJje3l6R3seCSt42MDBgeB6n98n7mWXLlmVt4+dBVAEAQPMCUQUAAKAh4VS9H/zgB/TII4/Q4YcfXtLP8ng8hjFHvKy2cWTM7nlalMxqm/l5AAAAmgvUVAEAAGg4du7cSeeffz7ddNNNtHr16lqfDgAAgCkOIlUAAAAaisHBQdE4l40p2BBi//79hv0ul4tmzpxZs/MDAAAw9YCoAgAA0FBwut+ePXvEMmvWrKz9bG3+9ttv1+TcAAAATE0cCnxeAQAATBHYfOKyyy4TSy0555xzhBX8Qw89VNPzAAAAUB5QUwUAAGBKwTbn7Oo3Ojpa9d/9zDPPiN/N9u8AAACaB0SqAAAATBk4ZZB7UjGHHnooOZ3VvbcYDofp/fffF49ZXMFmHQAAmgOIKgAAAAAAAAAoAaT/AQAAAAAAAEAJQFQBAAAAAAAAQAlAVAEAAAAAAABACUBUAQAAAAAAAEAJQFQBAAAAAAAAQAlAVAEAAAAAAABACUBUAQAAAAAAAEAJQFQBAAAAAAAAABXP/wcHxAm4+lSlwAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1000x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lens = SingletConfigurable(radius_of_curvature[0], refractive_index[0])\n",
    "lens.draw(num_rays=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will run the Monte Carlo simulation. In each of the 1000 iterations, we generate a random system, compute the rms spot radius and rms wavefront error, and record these values for later analysis.\n",
    "\n",
    "This code uses functionalities from the analysis and wavefront modules, which will be discussed in more detail in future tutorials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "rms_spot_radius = []\n",
    "rms_wavefront_error = []\n",
    "\n",
    "for k in range(num_systems):\n",
    "    # generate a random singlet\n",
    "    lens = SingletConfigurable(radius_of_curvature[k], refractive_index[k])\n",
    "\n",
    "    # get value of rms spot radius at field index = 0 and wavelength index = 0\n",
    "    spot = analysis.SpotDiagram(lens)\n",
    "    value = spot.rms_spot_radius()[0][0]\n",
    "    rms_spot_radius.append(value)\n",
    "\n",
    "    # get rms wavefront error\n",
    "    opd = wavefront.OPD(lens, field=(0, 0), wavelength=0.55)\n",
    "    rms_wavefront_error.append(opd.rms())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the output distributions of our computed metrics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGwCAYAAABPSaTdAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAANHBJREFUeJzt3QuYTeUex/H/jHEZBtOQmZExroVyKQqVS5JxOUWcU0qoHMkhhSTlFpVydDhJOs9zXFKJnKIIJaJkXJpSLiW3opoxjAZjMmbMOs//fZ69n9lzS2PP7L3f+X6eZ7Vnr7X22muv1ti/ea9BjuM4AgAAYKlgX58AAABAcSLsAAAAqxF2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYLcTXJ+APsrOz5ddff5XKlStLUFCQr08HAABcBB0q8MyZM1KzZk0JDi64/IawI2KCTkxMjK9PAwAAFMHRo0elVq1aBW4n7IiYEh3XxapSpYqvTwcAAFyE06dPm8IK1/d4QQg7Iu6qKw06hB0AAALLHzVBoYEyAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNVCfH0CtktKSpLU1FQJJOHh4RIVFeXr0wAAwCsIO8UcdLrfcaeknjkrgSS8ciVZ/cFyAg8AwAqEnWKkJToadC7vOEAqVouWQJCekijHNy4y507YAQDYgLBTAjTohEXG+vo0AAAolWigDAAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFbzadiZO3euNGvWTKpUqWKWtm3bypo1a9zbz507J8OGDZNq1apJWFiY9OnTR44dO+ZxjCNHjkiPHj2kYsWKUqNGDRkzZoxkZWX54NMAAAB/5NOwU6tWLXnhhRckISFBvvzyS+nUqZP07NlT9uzZY7aPHDlSVq5cKcuWLZNNmzbJr7/+Kr1793a//sKFCybonD9/XrZs2SKvv/66LFy4UCZOnOjDTwUAAPxJiC/f/Pbbb/d4/txzz5nSnq1bt5ogNG/ePFm8eLEJQWrBggXSuHFjs71Nmzby8ccfy969e+WTTz6RyMhIadGihUydOlXGjh0rkydPlnLlyvnokwEAAH/hN212tJRmyZIlcvbsWVOdpaU9mZmZ0rlzZ/c+jRo1ktq1a0t8fLx5ro9NmzY1QcclLi5OTp8+7S4dyk9GRobZJ+cCAADs5POws2vXLtMep3z58vLwww/L8uXLpUmTJpKUlGRKZsLDwz3212Cj25Q+5gw6ru2ubQWZNm2aVK1a1b3ExMQUy2cDAAC+5/Owc9VVV8nOnTtl27ZtMnToUBk4cKCpmipO48aNk1OnTrmXo0ePFuv7AQCAUtpmR2npTYMGDczPLVu2lB07dsi///1vufvuu03D49TUVI/SHe2NFRUVZX7Wx+3bt3scz9Vby7VPfrQUSRcAAGA/n5fs5JadnW3a1GjwKVu2rKxfv969bd++faarubbpUfqo1WDJycnufdatW2e6sWtVGAAAgE9LdrQ6qVu3bqbR8ZkzZ0zPq40bN8pHH31k2tIMGjRIRo0aJRERESbAPPLIIybgaE8s1aVLFxNq+vfvL9OnTzftdMaPH2/G5qHkBgAA+DzsaInMgAEDJDEx0YQbHWBQg85tt91mts+cOVOCg4PNYIJa2qM9rV599VX368uUKSOrVq0ybX00BFWqVMm0+ZkyZYoPPxUAAPAnPg07Oo5OYSpUqCBz5swxS0FiY2Nl9erVxXB2AADABn7XZgcAAMCbCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNcIOAACwmk/DzrRp0+T666+XypUrS40aNaRXr16yb98+j306duwoQUFBHsvDDz/ssc+RI0ekR48eUrFiRXOcMWPGSFZWVgl/GgAA4I9CfPnmmzZtkmHDhpnAo+Hkqaeeki5dusjevXulUqVK7v0GDx4sU6ZMcT/XUONy4cIFE3SioqJky5YtkpiYKAMGDJCyZcvK888/X+KfCQAA+Befhp21a9d6PF+4cKEpmUlISJD27dt7hBsNM/n5+OOPTTj65JNPJDIyUlq0aCFTp06VsWPHyuTJk6VcuXJ5XpORkWEWl9OnT3v1cwEAAP/hV212Tp06ZR4jIiI81r/11ltSvXp1ueaaa2TcuHGSnp7u3hYfHy9NmzY1QcclLi7OBJg9e/YUWH1WtWpV9xITE1NsnwkAAJTikp2csrOz5bHHHpObbrrJhBqXe++9V2JjY6VmzZry7bffmhIbbdfz3nvvme1JSUkeQUe5nuu2/GhgGjVqlPu5BiMCDwAAdvKbsKNtd3bv3i2bN2/2WP/QQw+5f9YSnOjoaLn11lvl4MGDUr9+/SK9V/ny5c0CAADs5xfVWMOHD5dVq1bJp59+KrVq1Sp039atW5vHAwcOmEdty3Ps2DGPfVzPC2rnAwAASg+fhh3HcUzQWb58uWzYsEHq1q37h6/ZuXOnedQSHtW2bVvZtWuXJCcnu/dZt26dVKlSRZo0aVKMZw8AAAJBiK+rrhYvXizvv/++GWvH1cZGGw2Hhoaaqird3r17d6lWrZppszNy5EjTU6tZs2ZmX+2qrqGmf//+Mn36dHOM8ePHm2NTVQUAAHxasjN37lzTA0sHDtSSGteydOlSs127jWuXcg00jRo1ktGjR0ufPn1k5cqV7mOUKVPGVIHpo5by3HfffWacnZzj8gAAgNIrxNfVWIXRHlI68OAf0d5aq1ev9uKZAQAAW/hFA2UAAIDiQtgBAABWI+wAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWK1IYeerr76SXbt2uZ+///770qtXL3nqqafk/Pnz3jw/AACAkg87Q4YMkR9++MH8fOjQIenbt69UrFhRli1bJk888cSlnREAAICvw44GnRYtWpifNeC0b99eFi9eLAsXLpR3333Xm+cHAABQ8mHHcRzJzs42P3/yySfSvXt383NMTIycOHHioo8zbdo0uf7666Vy5cpSo0YNUxW2b98+j33OnTsnw4YNk2rVqklYWJj06dNHjh075rHPkSNHpEePHqZ0SY8zZswYycrKKspHAwAAlilS2GnVqpU8++yz8sYbb8imTZtM0FCHDx+WyMjIiz6OvlaDzNatW2XdunWSmZkpXbp0kbNnz7r3GTlypKxcudKUIOn+v/76q/Tu3du9/cKFC+b9ta3Qli1b5PXXXzclTBMnTizKRwMAAJYJKcqLZs2aJf369ZMVK1bI008/LQ0aNDDr//e//8mNN9540cdZu3atx3MNKVoyk5CQYKrGTp06JfPmzTNVZJ06dTL7LFiwQBo3bmwCUps2beTjjz+WvXv3mhImDVpavTZ16lQZO3asTJ48WcqVK1eUjwgAAEpz2GnWrJlHbyyXf/7zn1KmTJkin4yGGxUREWEeNfRoaU/nzp3d+zRq1Ehq164t8fHxJuzoY9OmTT1KlOLi4mTo0KGyZ88eufbaa/O8T0ZGhllcTp8+XeRzBgAAlo6zk5qaKv/9739l3LhxcvLkSbNOS1iSk5OLdDxtA/TYY4/JTTfdJNdcc41Zl5SUZEpmwsPDPfbVYKPbXPvkrjpzPXftk19boapVq7oXbWsEAADsVKSw8+2330rDhg3lxRdflBkzZpjgo9577z0TfopC2+7s3r1blixZIsVNz1FLkVzL0aNHi/09AQBAAIWdUaNGyQMPPCD79++XChUquNdrr6zPPvvsTx9v+PDhsmrVKvn000+lVq1a7vVRUVGm4bErTLlobyzd5tond+8s13PXPrmVL19eqlSp4rEAAAA7FSns7NixwwwsmNsVV1xRYNVRQV3YNegsX75cNmzYIHXr1vXY3rJlSylbtqysX7/evU67pmtX87Zt25rn+qjth3JWn2nPLg0wTZo0KcrHAwAApb2BspaM5NeoVwcbvPzyy/9U1ZX2tNLpJnSsHVdQ0nY0oaGh5nHQoEGmJEkbLWuAeeSRR0zA0cbJSruqa6jp37+/TJ8+3Rxj/Pjx5th6ngAAoHQrUsnOHXfcIVOmTDE9pVRQUJApbdHu3jro38WaO3euaTPTsWNHiY6Odi9Lly517zNz5kz5y1/+Yo6r3dG1akrbBrlo7y+tAtNHDUH33XefDBgwwJwfAABAkUp2XnrpJfnrX/9qxsT5/fffpUOHDqZERcPGc88996eqsf6ItgmaM2eOWQoSGxsrq1evvuj3BQAApUeRwo5WL2m7mC+++EK++eYbSUtLk+uuu85jPBwErqysTDPBayDR4QkKapAOACjdihR2XHRMHF1gj4y0VElOSpIhI0abxuGBIrxyJVn9wXICDwDAO2FnxIgRZooIfczplVdekQMHDpjpJBCYss6lixNcRqp36C/h0bESCNJTEuX4xkVmiALCDgDAK2Hn3XfflQ8++CDPep0X64UXXiDsWCA0IkrCIgMj7AAA4PXeWCkpKabdTm7aNfzEiRPeOC8AAADfhR2twso9Y7las2aN1KtXzxvnBQAA4LtqLB3kT0c+Pn78uHTq1Mms01GOtUs6VVgAACDgw86DDz4oGRkZZkydqVOnmnV16tQxgwTqgH4AAAAB3/V86NChZtHSHZ3aISwszLtnBgAA4OtxdtSfmQsLAAAgIBooHzt2zEy8WbNmTQkJCTHzUuVcAAAAArpk5/777zcTf06YMMFM3KkTgQIAAFgTdjZv3iyff/65tGjRwvtnBAAA4OtqrJiYmIuasRwAACAgw46OpfPkk0/Kjz/+6P0zAgAA8HU11t133y3p6elSv359qVixYp7ZsU+ePOmt8wMAACj5sMMoyQAAwOqwM3DgQO+fCQAAgL+02VEHDx6U8ePHyz333CPJycnuiUD37NnjzfMDAAAo+bCzadMmadq0qWzbtk3ee+89SUtLM+u/+eYbmTRp0qWdEQAAgK/DjvbEevbZZ2XdunVSrlw593qdAX3r1q3ePD8AAICSDzu7du2SO++8M8/6GjVqyIkTJy7tjAAAAHwddsLDwyUxMTHP+q+//lquuOIKb5wXAACA78JO3759ZezYsZKUlGTmxcrOzpYvvvhCHn/8cRkwYIB3zgwAAMBXYef555+XRo0amWkjtHFykyZNpH379nLjjTeaHloAAAABO86OzomlJTovv/yyTJw40bTf0cBz7bXXSsOGDYvnLAEAAEoy7DRo0MCMp6PhRkt3AAAArKnGCg4ONiEnJSWleM4IAADA1212XnjhBRkzZozs3r3bm+cCAADgH3NjaY8rnfW8efPmZlDB0NBQj+3Meg4AAPwFs54DAACr/emwk5mZaebGmjBhgtStW7d4zgoAAMBXbXbKli0r7777rrfeHwAAwP8aKPfq1UtWrFjh/bMBAADwhzY72vV8ypQpZoqIli1bSqVKlTy2jxgxwlvnBwAAUPJhZ968eWYy0ISEBLPkpHNlEXYAAEBAh53Dhw97/0wAAAD8pc0OAACA1SU7Dz74YKHb58+fX9TzAQAA8H3Y+e233/KMvaNTR6SmpkqnTp28dW4AAAC+CTvLly/Psy47O1uGDh0q9evXv/SzAgAA8Lc2Ozob+qhRo2TmzJneOiQAAIB/NVA+ePCgZGVlefOQAAAAJV+NpSU4OTmOI4mJifLhhx/KwIEDL+2MAAAAfF2y8/XXX3ss3377rVn/0ksv/akZ0T/77DO5/fbbpWbNmmYwwtxTUNx///1mfc6la9euHvucPHlS+vXrJ1WqVDEDHQ4aNEjS0tKK8rEAAICFilSy8+mnn3rlzc+ePSvNmzc3Xdl79+6d7z4abhYsWOB+Xr58eY/tGnS0VGndunWmV9gDDzwgDz30kCxevNgr5wgAAErpCMraNkfnyMpp//79Zlb0OnXqXNRxunXrZpbCaLiJiorKd9t3330na9eulR07dkirVq3MutmzZ0v37t1lxowZpsQoPxkZGWZxOX369EWdLwAACDxFqsbS6qUtW7bkWb9t2zazzZs2btwoNWrUkKuuusp0bU9JSXFvi4+PN1VXrqCjOnfubHqG6bkUZNq0aVK1alX3EhMT49VzBgAAFrTZuemmm/Ksb9OmjezcuVO8RauwFi1aJOvXr5cXX3xRNm3aZEqCLly4YLYnJSWZIJRTSEiIREREmG0FGTdunJw6dcq9HD161GvnDAAALKjG0obCZ86cybNeg4MriHhD37593T83bdpUmjVrZgYt1NKeW2+9tcjH1aqx3G1/AACAnYpUstO+fXtTFZQz2OjPuu7mm2+W4lKvXj2pXr26HDhwwDzXtjzJycke+2hbIu2hVVA7HwAAULoUqWRHq5Q08Gg7mnbt2pl1n3/+uWnou2HDBikuP//8s2mzEx0dbZ63bdvWzMeVkJAgLVu2NOv0/XXqitatWxfbeQAAAMtLdpo0aWLG1rnrrrtMyYpWaQ0YMEC+//57ueaaay76ODoejrbxcbXz0V5e+vORI0fMtjFjxsjWrVvlxx9/NO12evbsKQ0aNJC4uDizf+PGjU27nsGDB8v27dvliy++kOHDh5vqr4J6YgEAgNKlSCU7SsPE888/f0lv/uWXX8ott9ySZ2RmHYV57ty5JlC9/vrrpvRG369Lly4ydepUj/Y2b731lgk42oZHe2H16dNHXn755Us6LwAAUMrDjg7yFxYWJn/729881i9btkzS09MvesqIjh07mqkmCvLRRx/94TG05xUDCAIAAK9WY2lDZG0onJt2A7/U0h4AAACfhx1tU1O3bt0862NjY802AACAgK7G0hIcbU+Te1qIb775RqpVq+atcwMuWlZWphw6dEgCiY7+zRAJAOCnYeeee+6RESNGSOXKlU0XdKWjGz/66KMeAwECJSEjLVWSk5JkyIjRZm62QBFeuZKs/mA5gQcA/DHsaI8o7Q6uPaB0egbXoILaMJk2OyhpWefSxQkuI9U79Jfw6FgJBOkpiXJ84yLT05CwAwB+GHbKlSsnS5culccff9yEntDQUDOdg7bZAXwlNCJKwiK5BwEAlxh29C/Rp59+2oSd3377zay77LLLTPXVs88+a9ohAAAABGTY0TmndIqGX375Rfr162dGMFZ79+6VhQsXmlGOt2zZYsIPAABAwIWdKVOmmCqsgwcPSmRkZJ5tOsKxPs6cOdPb5wkAAFD84+ysWLFCZsyYkSfoKG1kOX36dFm+fHnRzgQAAMDXYScxMVGuvvrqArfrJKBJSUneOC8AAICSDzs6RYT2viqIzlquc1UBAAAEZNiJi4szPbHOnz+fZ1tGRoZMmDBBunbt6s3zAwAAKNkGyq1atZKGDRvKsGHDpFGjRmbW8u+++05effVVE3jeeOONSzsjAAAAX4WdWrVqSXx8vPzjH/+QcePGmaCjgoKC5LbbbpNXXnlFYmJivHl+AAAAJTuooM52vmbNGjOg4P79+826Bg0a0FYHAADYM12E0oEDb7jhBu+eDQAAgC8bKAMAAAQawg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaiG+PgGgtMrKypRDhw5JoAkPD5eoqChfnwYAXDTCDuADGWmpkpyUJENGjJayZctKIAmvXElWf7CcwAMgYBB2AB/IOpcuTnAZqd6hv4RHx0qgSE9JlOMbF0lqaiphB0DAIOwAPhQaESVhkYETdgAgENFAGQAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaj4NO5999pncfvvtUrNmTQkKCpIVK1Z4bHccRyZOnCjR0dESGhoqnTt3lv3793vsc/LkSenXr59UqVLFjOw6aNAgSUtLK+FPAgAA/JVPw87Zs2elefPmMmfOnHy3T58+XV5++WV57bXXZNu2bVKpUiWJi4uTc+fOuffRoLNnzx5Zt26drFq1ygSohx56qAQ/BQAA8Gc+HVSwW7duZsmPlurMmjVLxo8fLz179jTrFi1aJJGRkaYEqG/fvvLdd9/J2rVrZceOHdKqVSuzz+zZs6V79+4yY8YMU2IEAABKN79ts3P48GFJSkoyVVcuVatWldatW0t8fLx5ro9adeUKOkr3Dw4ONiVBBcnIyJDTp097LAAAwE5+G3Y06CgtyclJn7u26WONGjU8toeEhEhERIR7n/xMmzbNBCfXEhMTUyyfAQAA+J7fhp3iNG7cODl16pR7OXr0qK9PCQAAlLaw45pR+dixYx7r9blrmz4mJyd7bM/KyjI9tAqbkbl8+fKm91bOBQAA2MlvZz2vW7euCSzr16+XFi1amHXatkbb4gwdOtQ8b9u2raSmpkpCQoK0bNnSrNuwYYNkZ2ebtj0AvC8rK1MOHTokgUTb9hX2BxAAu/k07Oh4OAcOHPBolLxz507T5qZ27dry2GOPybPPPisNGzY04WfChAmmh1WvXr3M/o0bN5auXbvK4MGDTff0zMxMGT58uOmpRU8swPsy0lIlOSlJhowYLWXLlpVAEV65kqz+YDmBByilfBp2vvzyS7nlllvcz0eNGmUeBw4cKAsXLpQnnnjCjMWj4+ZoCc7NN99suppXqFDB/Zq33nrLBJxbb73V9MLq06ePGZsHgPdlnUsXJ7iMVO/QX8KjYyUQpKckyvGNi8y/IYQdoHTyadjp2LGjGU+nIDqq8pQpU8xSEC0FWrx4cTGdIYD8hEZESVhkYIQdAPDbBsoAAADeQNgBAABWI+wAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYLUQX58AABS3rKxMOXTokASS8PBwiYqK8vVpAFYg7ACwWkZaqiQnJcmQEaOlbNmyEijCK1eS1R8sJ/AAXkDYAWC1rHPp4gSXkeod+kt4dKwEgvSURDm+cZGkpqYSdgAvIOwAKBVCI6IkLDIwwg4A76KBMgAAsBphBwAAWM2vw87kyZMlKCjIY2nUqJF7+7lz52TYsGFSrVo1CQsLkz59+sixY8d8es4AAMC/+HXYUVdffbUkJia6l82bN7u3jRw5UlauXCnLli2TTZs2ya+//iq9e/f26fkCAAD/4vcNlENCQvLtjXDq1CmZN2+eLF68WDp16mTWLViwQBo3bixbt26VNm3aFHjMjIwMs7icPn26mM4eAAD4mt+X7Ozfv19q1qwp9erVk379+smRI0fM+oSEBMnMzJTOnTu799Uqrtq1a0t8fHyhx5w2bZpUrVrVvcTExBT75wAAAL7h12GndevWsnDhQlm7dq3MnTtXDh8+LO3atZMzZ85IUlKSlCtXzowymlNkZKTZVphx48aZkiHXcvTo0WL+JAAAwFf8uhqrW7du7p+bNWtmwk9sbKy88847EhoaWuTjli9f3iwAAMB+fl2yk5uW4lx55ZVy4MAB047n/PnzZoTRnLQ3FiOOAgCAgAw7aWlpcvDgQYmOjpaWLVuaeW7Wr1/v3r5v3z7Tpqdt27Y+PU8AAOA//Loa6/HHH5fbb7/dVF1pt/JJkyZJmTJl5J577jENiwcNGiSjRo2SiIgIqVKlijzyyCMm6BTWEwsAAJQufh12fv75ZxNsUlJS5PLLL5ebb77ZdCvXn9XMmTMlODjYDCaoXcnj4uLk1Vdf9fVpAwAAP+LXYWfJkiWFbq9QoYLMmTPHLAAAAAHfZgcAAODPIuwAAACrEXYAAIDVCDsAAMBqft1AGQBKq6ysTDl06JAE2sCvDOoKf0TYAQA/k5GWKslJSTJkxGgzeGqgCK9cSVZ/sJzAA79D2AEAP5N1Ll2c4DJSvUN/CY+OlUCQnpIoxzcuMlP4EHbgbwg7AOCnQiOiJCwyMMIO4M9ooAwAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsFuLrEwAA2CErK1MOHTokgSQ8PFyioqJ8fRooZoQdAMAly0hLleSkJBkyYrSULVtWAkV45Uqy+oPlBB7LEXYAAJcs61y6OMFlpHqH/hIeHSuBID0lUY5vXCSpqamEHcsRdgAAXhMaESVhkYERdlB60EAZAABYjbADAACsRtgBAABWo80OAKDUort86UDYAQCUSnSXLz0IOwCAUonu8qUHYQcAUKrRXd5+hB0AAAII7Yz+PMIOAAABgnZGRUPYAQAgQNDOqGgIOwAABBjaGf05DCoIAACsRtgBAABWI+wAAACrWRN25syZI3Xq1JEKFSpI69atZfv27b4+JQAA4AesCDtLly6VUaNGyaRJk+Srr76S5s2bS1xcnCQnJ/v61AAAgI9ZEXb+9a9/yeDBg+WBBx6QJk2ayGuvvSYVK1aU+fPn+/rUAACAjwV81/Pz589LQkKCjBs3zr0uODhYOnfuLPHx8fm+JiMjwywup06dMo+nT5/26rmlpaVJ9oULcibxsGRlpEsgOHv8Z3GcbElL+lHKyAUJBJxzyQnE8+acSwbnXDIC8Zx/P3nMfBfqd6K3v2ddx3Mcp/AdnQD3yy+/6Cd0tmzZ4rF+zJgxzg033JDvayZNmmRew8LCwsLCwiIBvxw9erTQrBDwJTtFoaVA2sbHJTs7W06ePCnVqlWToKAgsZUm4JiYGDl69KhUqVLF16fjF7gmeXFN8uKa5MU1yYtrUvLXQ0t0zpw5IzVr1ix0v4APO9WrV5cyZcrIsWPHPNbr84KGpS5fvrxZck9SVlroTccvoieuSV5ck7y4JnlxTfLimpTs9ahatar9DZTLlSsnLVu2lPXr13uU1Ojztm3b+vTcAACA7wV8yY7SKqmBAwdKq1at5IYbbpBZs2bJ2bNnTe8sAABQulkRdu6++245fvy4TJw4UZKSkqRFixaydu1aiYyM9PWp+RWtutOxiHJX4ZVmXJO8uCZ5cU3y4prkxTXx3+sRpK2UfX0SAAAAxSXg2+wAAAAUhrADAACsRtgBAABWI+wAAACrEXYCyJw5c6ROnTpSoUIFad26tWzfvr3Q/ZctWyaNGjUy+zdt2lRWr17t3paZmSljx4416ytVqmRGnxwwYID8+uuvHsfQ99NRpXMuL7zwgth4TdT999+f5/N27drVYx8dbbtfv35mkCwdjHLQoEFmzhdbr0nu6+Fa/vnPf1p5n+zZs0f69Onj/kw6lEVRjnnu3DkZNmyYGZk9LCzMHDP34Kc2XZNp06bJ9ddfL5UrV5YaNWpIr169ZN++fR77dOzYMc998vDDD4ut12Ty5Ml5Pq/+rgXKfTLHy9cjv38ndNHPX+z3iDfnqULxWbJkiVOuXDln/vz5zp49e5zBgwc74eHhzrFjx/Ld/4svvnDKlCnjTJ8+3dm7d68zfvx4p2zZss6uXbvM9tTUVKdz587O0qVLne+//96Jj483c4m1bNnS4zixsbHOlClTnMTERPeSlpbm2HhN1MCBA52uXbt6fN6TJ096HEe3N2/e3Nm6davz+eefOw0aNHDuuecex9ZrkvNa6KLHDgoKcg4ePGjlfbJ9+3bn8ccfd95++20nKirKmTlzZpGO+fDDDzsxMTHO+vXrnS+//NJp06aNc+ONNzq2XpO4uDhnwYIFzu7du52dO3c63bt3d2rXru1xH3To0MG8V8775NSpU46t10TnYbz66qs9Pu/x48c99vHX+2RJMVyP5ORkj2uxbt06M6/Vp59+Wuz3CGEnQGgQGTZsmPv5hQsXnJo1azrTpk3Ld/+77rrL6dGjh8e61q1bO0OGDCnwPfRm1Rvvp59+8vgSy++mtfWaaNjp2bNnge+pgUCv0Y4dO9zr1qxZY778dVLa0nCf6PXp1KmTxzqb7pOL+Vx/dEz9Y0JD47Jly9z7fPfdd+be0T8sbLwm+X2x6efdtGmTxxfZo48+6vij4rgmGnb0D6OC+PN9ckMJ3CN6L9SvX9/Jzs4u9nuEaqwAcP78eUlISJDOnTu71wUHB5vn8fHx+b5G1+fcX8XFxRW4vzp16pQpMsw9T5hWR2gR67XXXmuqLrKyssTma7Jx40ZTDH/VVVfJ0KFDJSUlxeMYen10tG4XPaa+97Zt28T2+0SL1z/88ENTdZebLfeJN46p27WqOOc+Wn1Ru3btIr+vP1+Tgv49URERER7r33rrLTOn4TXXXGMmZU5PTxdfK85rsn//ftNMoF69eqb6+8iRI+5t/nqfnC+Be0Tf480335QHH3wwzwTcxXGPWDGCsu1OnDghFy5cyDMitD7//vvv832NjiSd3/66Pj9ab6xteO655x6PCdtGjBgh1113nfkHa8uWLebGS0xMlH/9619i4zXR9jm9e/eWunXrysGDB+Wpp56Sbt26mV9wnXBW99UglFNISIi5PgVdW5vuk9dff920ydBrlJNN94k3jqnXT+fty/2HQ2HXNpCvSW46P+Fjjz0mN910k/nCcrn33nslNjbWfPl/++235t8cbdfz3nvviY3XRNu5LFy40PzhpL8PzzzzjLRr1052795tfo/89T45UQL3yIoVKyQ1NdW0k8ypuO4Rwg7MXxZ33XWXVmnK3Llz88w75tKsWTPzizlkyBDTGNEfhgD3tr59+7p/1sa6+pnr169vSntuvfVWKe3mz59v/jrVBoul+T5B4bTBqX6hb9682WP9Qw895PH7FR0dbX6v9A8L/T2zjf6hlPP3QsOPfpG/8847+ZaOlibz5s0z10dDTUncI1RjBQAtztNShdwt9PV5VFRUvq/R9Rezvyvo/PTTT7Ju3TqPUp386C+rVk/8+OOPYus1yUmLnvW9Dhw44D5GcnKyxz56PbSHVmHHseGafP755+YvrL///e9/eC6BfJ9445j6qMX0+pert97Xn69JTsOHD5dVq1bJp59+KrVq1frD+0S5fr9svSYuWoJz5ZVXevx74o/3SfVivh76ffPJJ59c9L8l3rhHCDsBQP9Kbtmypaxfv96jmFift23bNt/X6Pqc+ysNMzn3dwUdrVPWG0/bW/yRnTt3mrrb3FU5tlyT3H7++WfTZkf/unAdQ/9h0vpslw0bNpj3dv1S2npN9C8xPX7z5s2tvk+8cUzdXrZsWY99NChqe42ivq8/XxOlJcMadJYvX25+J7Qq+GLuE+X6/bLtmuSmQ1RoCYXr8/rrfVKumK/HggULzL8NPXr0KLl7xOtNnlFs3QDLly/vLFy40PQIeuihh0w3wKSkJLO9f//+zpNPPunRpTgkJMSZMWOGad2vvQJydik+f/68c8cddzi1atUy3URzdvPLyMgw+2zZssW0qNft2s34zTffdC6//HJnwIABjo3X5MyZM6brpPaCOHz4sPPJJ5841113ndOwYUPn3LlzHl3Pr732Wmfbtm3O5s2bzXZ/6nruzWviol0/K1as6MydOzfPe9p2n+j9//XXX5slOjra3BP68/79+y/6mK4uxdr1esOGDaZLcdu2bc1i6zUZOnSoU7VqVWfjxo0e/56kp6eb7QcOHDDDE+i10N+v999/36lXr57Tvn17x9ZrMnr0aHM99PPq75oO91G9enXTU83f75MlxXA9XL269POOHTs2z3sW5z1C2Akgs2fPNjeJjn2g3QJ1nJec3fW023RO77zzjnPllVea/XWshw8//NC9TW8kzbr5La4xDxISEkw3ZP0HrEKFCk7jxo2d559/3uOL36Zrov8od+nSxXxR6xe+dp/U8R5yfoGplJQUE27CwsKcKlWqOA888IAJSjZeE5f//Oc/TmhoqOkqm5tt90lBvxu638UeU/3+++/OP/7xD+eyyy4zQfHOO+80X/62XpOC/j3RsXfUkSNHzJdWRESE+RLV8anGjBnjN+PsFMc1ufvuu80Xvx7viiuuMM/1Cz1Q7pPZxfB789FHH5n1+/bty/N+xXmPBOl/Lq1sCAAAwH/RZgcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAKERQUJCtWrDA/68Sm+tw1X09x0TmIGjduLBcuXJCS9tprr8ntt99e4u8LFCfCDhAg7r//fvNFq4tOHqgTLT7xxBNy7tw5j/1c+2zdutVjfUZGhpnsVbdt3LjRvX7Tpk3SqVMniYiIkIoVK0rDhg1l4MCBZjbmgnzzzTdyxx13mMn8KlSoIHXq1JG77747z4zwl0qPO2vWrIvaz/W59TM0bdpU/vvf/4q3xcTESGJiolxzzTVSnPT/6/jx483M0yXtwQcflK+++srMcg/YgrADBJCuXbuaL9tDhw7JzJkz5T//+Y9MmjQp3y9lnVk4J52NOiwszGPd3r17zTFbtWoln332mezatUtmz55tZj0uqFTh+PHjcuutt5pw9NFHH8l3331n3qtmzZpy9uxZ8ZUpU6aYa7N792657777ZPDgwbJmzRqvvoeGj6ioKAkJCZHisnnzZjMzdp8+fcQX9P/9vffeKy+//LJP3h8oFpc8uxaAEqGT7vXs2dNjXe/evc0M7Dnpr/X48ePNJKWuGafVbbfd5kyYMMFjsledrbxOnTp/6jyWL19uZkrPzMwscB89vr7PqlWrnKZNm5pJ/XSy0Nyzqf/vf/9zmjRpYiYa1IlXdfZ1F51AMPekggXR1+pnyUknExw5cqT7+fbt282s09WqVTPXRicc1ElMc/rhhx+cdu3amfPVCU0//vhj8776mXNOdqizOSud5FInQM19fXKeq84G37FjRzNxbOXKlZ3rrrvO2bFjR4GfZdiwYc5f//pXj3U6G33z5s2defPmOTExMU6lSpXMLONZWVnOiy++6ERGRpoJbJ999lmP1+l5vPbaa06PHj3MRK6NGjUys9TrTNR6fXXiSZ1hO/fklJs2bTL/T3LeP0Ago2QHCFBagrFlyxbzl3huLVu2NFU77777rnl+5MgRU3LTv39/j/20lEJLQ3TbxdLXZGVlmZKiP5pHeMyYMfLSSy/Jjh075PLLLzdtQTIzM822hIQEueuuu6Rv376mRGny5MkyYcIEWbhwodn+3nvvSa1atdwlNrpcjOzsbPO5f/vtN49rc+bMGVM9pyUnWsWn1XXdu3c3612v6927t3nNtm3bTNuVsWPHyqXq16+f+Rx6DfQzP/nkk6YasiBafaQlbblpaY+WVK1du1befvttmTdvnvTo0UN+/vlnUxX54osvmqovPfecpk6dKgMGDDDtjBo1amRKbYYMGSLjxo2TL7/80vw/HD58uMdr9P31/3HuYwEBy9dpC8DFl+yUKVPG/FWvJQ/66xscHGxKR3JylUTMmjXLueWWW8y6Z555xrnzzjud3377zaNkR0sG7r//frMuKirK6dWrlzN79mzn1KlThZ7LU089ZUp3tPSka9euzvTp052kpKQ8JTtLlixxr0tJSTGlC0uXLjXP7733XlPalNOYMWNMSU9hJTb50f20JEKvjZ6Xvreem5ZgFOTChQumpGXlypXm+UcffWRe+8svv7j3WbNmzSWX7Oh7LFy40LlYerxFixblKdnRUpjTp0+718XFxZlSOf0cLldddZUzbdq0PKV8LvHx8WadlhC5vP32206FChXynMdll132p84b8GeU7AAB5JZbbjF/oetf3FpK8cADDxTYtkPbrcTHx5v2PVpaog1P82uDou1ttHRg+vTpcsUVV8jzzz8vV199daElKc8995wkJSWZ0g/dVx+11EBLaHJq27at+2dt43PVVVeZNj5KH2+66SaP/fX5/v37i9QLSUuR9Nps2LBBWrdubdo0NWjQwL392LFjph2PluhUrVpVqlSpImlpaabUy3U+2tZJ2x7ld/5FNWrUKPn73/8unTt3lhdeeMGU0BTm999/N42+c9OSusqVK7ufR0ZGSpMmTSQ4ONhjXe5G4s2aNfPYrrQBd8512sj99OnTHq8LDQ2V9PT0P/VZAX9F2AECSKVKlcwXePPmzWX+/Pkm9Gh1Rn6059Vf/vIXGTRokPky69atW4HH1ZCjVVyvvPKK7Nmzx+yvAaYwevy//e1vMmPGDBMUNCToz75SvXp1c23atWsny5YtkxEjRpgG2C4aDjUM/fvf/zbVf/qzfobCep39EQ0auavyXNV0Llo9p9dUq5w0iGlA0SrAwj6HVsHllrvqy9UrL/c6rY4r6HW6vaB1uV938uRJU/UI2ICwAwQo/aJ96qmnTDsNLQ3Ij5bmaDdzbbNxsd2YL7vsMomOjv5TPau0nUv9+vXzvCZn93f9Av/hhx/M+DFKH7/44guP/fX5lVde6T7XwnqFFUZLaLQrvLZLyXlsDUDaTkdLo8qXLy8nTpxwb9fzOXr0qEeJVu7u+7lpGNA2Pzk/d35j8OhnGjlypHz88cemXVDunnI5XXvttR4hzRe09EkDr54LYAPCDhDAtGRFg8GcOXPy3a7dyrWruDbyzY92XR86dKj5EtYvOC2B0Ea5+ljQwHKrVq0yVWT6qOFl3759pkRn9erV0rNnT4999X11gDxtTK3jBGmpRa9evcy20aNHm23agFaP8/rrr5uSpccff9yj6kYbT//yyy8eweRiPProo7Jy5UrTCFdp9dUbb7xhSqG0REwbDmtVjYtWM2ko0RIgHUdIGwo//fTThb6HVpfpuD4aOvX6LV682N3AWmkI1ca/Gjh/+uknE7i0obIr8OUnLi7ONKL2Jf3s9erVMwEWsAFhBwhgOt6Lfplqe5v8SmK0ikIDRn49ttQNN9xg2q08/PDDprSjQ4cOpjRDRwzWn/Oj1TD6Ba9hpUWLFtKmTRt55513zCB+uXt7aRsVDR3aO0zb+Gj4cJ3LddddZ163ZMkSM0jfxIkTTTjSUOSiz3XUYv3S/bNVKnqeXbp0McdVWt2npUv6vnqeWsqjgyLmLCnT6iUNKHpdtJ2Ntk0qjLZDevPNN03Q03Yw2ktKq61cNIimpKSYkjUNUtr7TKsTn3nmmQKPqSFMw6aGSF/Rz6HtmwBbBGkrZV+fBAC7aEmGNqbWcBEeHu7r0wk42thaGwxryVtJ06ClI2praZs25AZsQMkOAPgZrT6LjY3N02i4JGibpUWLFhF0YBVKdgB4HSU7APwJYQcAAFiNaiwAAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAQGz2f616yMqLZb5KAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(rms_spot_radius, color=\"C0\", edgecolor=\"k\", alpha=0.8)\n",
    "plt.xlabel(\"RMS Spot Radius (mm)\")\n",
    "plt.ylabel(\"Occurrences\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGwCAYAAABPSaTdAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAANLNJREFUeJzt3QucTfX+//GPMWPGbYxLjMk11FAuJUkk4biVEv9TpBDl5FCkpBFdcChdOErUSdQ5OYly6UIxREKJlEvJZUQ1Y9zGjDTMZf0en+/57/2YPTeMPfbe3/16Ph7L2Guvvfd37bVm9nt/L+tbwnEcRwAAACwV4usCAAAAFCfCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1UJ9XQB/kJ2dLb///ruUL19eSpQo4eviAACAc6CXCkxLS5OYmBgJCSm4/oawI2KCTs2aNX1dDAAAUAQHDx6UGjVqFHg/YUfE1Oi43qzIyEhfFwcAAJyD1NRUU1nh+hwvCGFHxN10pUGHsAMAQGA5WxcUOigDAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWC3U1wWwXVJSkqSkpEggiYqKkujoaF8XAwAAryDsFHPQ+X+3dZX0k8clkESUqygLly4j8AAArEDYKUZao6NBZ0KH8lK3SoQEgoQj6TIu/rgpO2EHAGADws5FoEEnNrqMBI40XxcAAACvoYMyAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1XwadiZPniwtWrSQ8uXLS9WqVaVHjx6ya9cuj23atWsnJUqU8FgefPBBj20OHDggt9xyi5QpU8Y8z6hRoyQzM/Mi7w0AAPBHob588TVr1sjQoUNN4NFwMmbMGOnUqZPs3LlTypYt697ugQcekPHjx7tva6hxycrKMkEnOjpa1q9fL4mJidKvXz8JCwuTSZMmXfR9AgAA/sWnYWf58uUet+fOnWtqZjZv3ixt27b1CDcaZvLz+eefm3C0cuVKqVatmjRr1kwmTJggo0ePlmeeeUZKlSqV5zGnT582i0tqaqpX9wsAAPgPv+qzc+LECfOzUqVKHuvfffddqVKlilx11VUSFxcnp06dct+3YcMGady4sQk6Lp07dzYBZseOHQU2n1WoUMG91KxZs9j2CQAABHHNTk7Z2dkyYsQIad26tQk1LnfffbfUrl1bYmJi5IcffjA1Ntqv58MPPzT3JyUleQQd5bqt9+VHA9PIkSPdtzUYEXgAALCT34Qd7buzfft2Wbduncf6wYMHu/+vNTjVq1eXDh06yN69e6VevXpFeq3w8HCzAAAA+/lFM9awYcPk448/ltWrV0uNGjUK3bZly5bm5549e8xP7ctz6NAhj21ctwvq5wMAAIKHT8OO4zgm6CxatEhWrVoldevWPetjtm7dan5qDY9q1aqVbNu2TZKTk93brFixQiIjI6VRo0bFWHoAABAIQn3ddDVv3jxZsmSJudaOq4+NdhouXbq0aarS+7t16yaVK1c2fXYeeeQRM1KrSZMmZlsdqq6h5t5775UpU6aY5xg7dqx5bpqqAACAT2t2Zs6caUZg6YUDtabGtcyfP9/cr8PGdUi5BprY2Fh59NFHpVevXvLRRx+5n6NkyZKmCUx/ai3PPffcY66zk/O6PAAAIHiF+roZqzA6QkovPHg2Olrr008/9WLJAACALfyigzIAAEBxIewAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVCDsAAMBqob4uAPxPRmaW7Nu3TwJJVFSUREdH+7oYAAA/5NOwM3nyZPnwww/lp59+ktKlS8sNN9wgzz//vFxxxRXubdLT0+XRRx+V9957T06fPi2dO3eW1157TapVq+be5sCBAzJkyBBZvXq1lCtXTvr372+eOzSULHe+DqdlSGJSssQNHyxhYYHz/kWUqygLly4j8AAA8vDpp9maNWtk6NCh0qJFC8nMzJQxY8ZIp06dZOfOnVK2bFmzzSOPPCKffPKJLFiwQCpUqCDDhg2Tnj17yldffWXuz8rKkltuucV8yK1fv14SExOlX79+EhYWJpMmTfLl7gWktPRMCQvJlmfbl5XYmEgJBAlH0mVc/HFJSUkh7AAA/CvsLF++3OP23LlzpWrVqrJ582Zp27atnDhxQmbPni3z5s2T9u3bm23mzJkjDRs2lI0bN8r1118vn3/+uQlHK1euNLU9zZo1kwkTJsjo0aPlmWeekVKlSvlo7wJbncrhEhtdRgJHmq8LAADwU37VQVnDjapUqZL5qaEnIyNDOnbs6N4mNjZWatWqJRs2bDC39Wfjxo09mrW0qSs1NVV27NiR7+toc5jen3MBAAB28puwk52dLSNGjJDWrVvLVVddZdYlJSWZmhntfJqTBhu9z7VNzqDjut91X360P482ibmWmjVrFtNeAQAAX/ObsKN9d7Zv3246Ihe3uLg4U4vkWg4ePFjsrwkAAHzDL4bbaKfjjz/+WNauXSs1atRwr9fOpmfOnDEdT3PW7hw6dMjdEVV/fvPNNx7Pp/e77stPeHi4WQAAgP18WrPjOI4JOosWLZJVq1ZJ3bp1Pe5v3ry5GVUVHx/vXrdr1y4z1LxVq1bmtv7ctm2bJCcnu7dZsWKFREZGSqNGjS7i3gAAAH8U6uumKx1ptWTJEilfvry7j432o9Hr7ujPQYMGyciRI02nZQ0wDz30kAk4OhJL6VB1DTX33nuvTJkyxTzH2LFjzXNTewMAAHwadmbOnGl+tmvXzmO9Di8fMGCA+f/UqVMlJCREevXq5XFRQZeSJUuaJjC9qKCGIL0+j15UcPz48Rd5bwAAgD8K9XUz1tlERETIjBkzzFKQ2rVry6effurl0gEAABv4zWgsAACA4kDYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxUp7GzZskW2bdvmvr1kyRLp0aOHjBkzRs6cOePN8gEAAFyQ0KI86G9/+5s88cQT0rhxY9m3b5/07t1b7rjjDlmwYIGcOnVKpk2bdmGlAs5TRmaWORcDSVRUlERHR/u6GABgvSKFnZ9//lmaNWtm/q8Bp23btjJv3jz56quvTPAh7OBiOpyWIYlJyRI3fLCEhRXplPaJiHIVZeHSZQQeAChmRfpkcBxHsrOzzf9Xrlwpt956q/l/zZo15ciRI94tIXAWaemZEhaSLc+2LyuxMZESCBKOpMu4+OOSkpJC2AEAfww71157rUycOFE6duwoa9askZkzZ5r1CQkJUq1aNW+XETgndSqHS2x0GQkcab4uAAAEhSJ1UNZmKu2kPGzYMHnyySelfv36Zv3ChQvlhhtu8HYZAQAALm7NTpMmTTxGY7m88MILUrJkyaKXBgAAwF+us6N9Dd58802Ji4uTY8eOmXU7d+6U5ORkb5YPAADg4tfs/PDDD9KhQwczdHb//v3ywAMPSKVKleTDDz+UAwcOyDvvvHNhpQIAAPBlzc7IkSPlvvvuk927d0tERIR7fbdu3WTt2rXeKhsAAIBvws6mTZvMhQVzu/TSSyUpKenCSwUAAODLsBMeHi6pqan5Xmzwkksu8Ua5AAAAfBd2brvtNhk/frxkZGSY2yVKlDB9dUaPHi29evXyTskAAAB8FXZeeuklOXnypFStWlX+/PNPuemmm8y1dsqXLy//+Mc/vFEuAAAA343GqlChgqxYscLMhfX999+b4HPNNdeYKyoDAAD4kwuaNbF169ZmAQAAsKoZ6+GHH5bp06fnWf/qq6/KiBEjvFEuAAAA34WdDz74IN8aHZ0XS+fHAgAA8BdFCjtHjx41/XZyi4yMlCNHjnijXAAAAL4LOzryavny5XnWL1u2TC677DJvlAsAAMB3HZR1uohhw4bJ4cOHpX379mZdfHy8GZI+bdo075QMAADAV2Fn4MCBcvr0aXNNnQkTJph1derUkZkzZ0q/fv28US4AAADfDj0fMmSIWbR2p3Tp0lKuXDnvlAgAAMBfrrOjmAsLAABY10H50KFDcu+990pMTIyEhoZKyZIlPRYAAICArtkZMGCAmfhz3LhxUr16dTMRKAAAgDVhZ926dfLll19Ks2bNvF8iAAAAXzdj1axZUxzHueAXX7t2rXTv3t00h2nt0OLFi/PUIOn6nEuXLl08tjl27Jj07dvXXNAwKipKBg0aZCYmBQAAKHLY0WvpPPHEE7J///4Lehf/+OMPadq0qcyYMaPAbTTcJCYmupf//ve/Hvdr0NmxY4eZhf3jjz82AWrw4MEXVC4AABDkzVh33XWXnDp1SurVqydlypSRsLCwPLUt56Jr165mKUx4eLhER0fne9+PP/5oruS8adMmufbaa826V155Rbp16yYvvviiqTECAADBrUhh52JeJfmLL76QqlWrSsWKFc3VmidOnCiVK1c2923YsME0XbmCjurYsaOEhITI119/LXfccUe+z6kXRNTFJTU19SLsCQAACJiw079/f7kYtAmrZ8+eUrduXdm7d6+MGTPG1ARpyNEh7klJSSYI5aRD4StVqmTuK8jkyZPl2WefvQh7AAAAArLPjtLwMXbsWOnTp48kJye7JwLV/jPe0rt3b7ntttukcePG0qNHD9MnR5ustLbnQsTFxcmJEyfcy8GDB71WZgAAYEHYWbNmjQkg2lT04Ycfukc/ff/99/L0009LcdEZ1atUqSJ79uwxt7UvjytouWRmZpo+QwX183H1A9LRWzkXAABgpyKFHR2JpX1ndARUqVKl3Ou1T83GjRuluPz6669y9OhRcyFD1apVK0lJSZHNmze7t1m1apVkZ2dLy5Yti60cAADA8j4727Ztk3nz5uVZr/1njhw5cs7PozVCrloalZCQIFu3bjV9bnTRfjW9evUytTTabPb4449L/fr1pXPnzmb7hg0bmn49DzzwgMyaNUsyMjJk2LBhpvmLkVgAAKDINTs6AkqveZPbd999J5deeuk5P8+3334rV199tVnUyJEjzf+feuop0wH5hx9+MH12Lr/8cnOxwObNm5srN2szlMu7774rsbGx0qFDBzPkvE2bNvLGG29wdAEAQNFrdrTmZPTo0bJgwQJzVWNtNvrqq6/ksccek379+p3z87Rr167QKzF/9tlnZ30OrQHKr5YJAACgyDU7kyZNMrUpOm2ENkU1atRI2rZtKzfccIMZoQUAABCwNTtaE6PXsJk+fbppbtL+Oxp4tPmpQYMGxVNKAACAixl2tJOwXk9Hw43W7gAAAFjTjKVTMWjI0SHgAAAAVvbZee6552TUqFGyfft275cIAADA16OxdMSVznretGlTc1HB0qVLF2nWcwAAAAn2Wc8BAAAuatjRqxTr3Fjjxo0zs5EDAABY1WcnLCxMPvjgg+IpDQAAgD90UO7Ro4csXrzY22UBAADwjz47OvR8/PjxZooIna+qbNmyHvc//PDD3iofAADAxQ87s2fPNpOBbt682Sw56VxZhB0AABDQYSchIcH7JQEAAPCXPjsAAABW1+wMHDiw0PvfeuutopYHAADA92Hn+PHjea69o1NHpKSkSPv27b1VNgAAAN+EnUWLFuVZl52dLUOGDJF69epdeKkAAAD8rc+OzoY+cuRImTp1qreeEgAAwL86KO/du1cyMzO9+ZQAAAAXvxlLa3BychxHEhMT5ZNPPpH+/ftfWIkAAAB8HXa+++67PE1Yl1xyibz00ktnHakFAADg92Fn9erV3i8JAACAv/TZ0Sso7969O896Xbd//35vlAsAAMArihR2BgwYIOvXr8+z/uuvvzb3AQAABHTY0T47rVu3zrP++uuvl61bt3qjXAAAAL4LOzqzeVpaWp71J06ckKysLG+UCwAAwHdhp23btjJ58mSPYKP/13Vt2rTxTskAAAB8NRrr+eefN4HniiuukBtvvNGs+/LLLyU1NVVWrVrljXIBAAD4rmanUaNG8sMPP8idd94pycnJpkmrX79+8tNPP8lVV13lnZIBAAD4qmZHxcTEyKRJk7xRBgAAAP+q2ZkzZ44sWLAgz3pd9/bbb3ujXAAAAL4LO9oRuUqVKnnWV61aldoeAAAQ+GHnwIEDUrdu3Tzra9eube4DAAAI6LCjNTjaQTm377//XipXruyNcgEAAPgu7PTp00cefvhhMyGoXl9HFx1yPnz4cOndu7d3SgYAAOCr0VgTJkwwE3526NBBQkP/9xQaePr370+fHQAAEPhhp1SpUjJ//nx57LHHTOgpXbq0NG7c2PTZAQAACOiwk5KSIk8++aQJO8ePHzfrKlasaJqvJk6cKFFRUcVRTgAAgOIPO8eOHZNWrVrJb7/9Jn379pWGDRua9Tt37pS5c+dKfHy8rF+/3oQfAACAgAs748ePN01Ye/fulWrVquW5r1OnTubn1KlTvV1OAACA4h+NtXjxYnnxxRfzBB0VHR0tU6ZMkUWLFhWtJAAAAL6u2UlMTJQrr7yywPt1EtCkpCRvlAuwXkZmluzbt08CjfbL0y83AGBl2NEpInT0VY0aNfK9PyEhQSpVquStsgHWOpyWIYlJyRI3fLCEhRV5Pl6fiChXURYuXUbgARAwzuuvbOfOnc1IrBUrVpi+OzmdPn1axo0bJ126dPF2GQHrpKVnSlhItjzbvqzExkRKoEg4ki7j4o+bUZmEHQDWdlC+9tprpUGDBjJ06FCJjY0Vx3Hkxx9/lNdee80Enn//+9/FV1rAMnUqh0tsdBkJLGm+LgAAFF/Y0earDRs2yN///neJi4szQUeVKFFC/vKXv8irr74qNWvWPL8SAAAAFKPz7iygs50vW7bMXFBw9+7dZl39+vXpqwMAAPxSkXtG6oUDr7vuOu+WBgAAwB9mPQcAAAgUhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFbzadhZu3atdO/eXWJiYsyFCXVW9Zz0ooVPPfWUVK9eXUqXLi0dO3Z0X9vH5dixY9K3b1+JjIw0ExQOGjRITp48eZH3BAAA+Cufhp0//vhDmjZtKjNmzMj3/ilTpsj06dNl1qxZ8vXXX0vZsmXN/Fzp6enubTTo7Nixw8zX9fHHH5sANXjw4Iu4FwAAwJ/5dLrlrl27miU/Wqszbdo0GTt2rNx+++1m3TvvvCPVqlUzNUC9e/c2c3ItX75cNm3aZObsUq+88op069ZNXnzxRVNjlB+dw0sXl9TU1GLZPwAA4Ht+22cnISFBkpKSTNOVS4UKFaRly5Zmfi6lP7XpyhV0lG4fEhJiaoIKMnnyZPNcroX5vAAAsJffhh0NOkprcnLS26779GfVqlU97g8NDTXzdLm2yY9OYnrixAn3cvDgwWLZBwAAEOTNWL4SHh5uFgAAYD+/rdmJjo42Pw8dOuSxXm+77tOfycnJHvdnZmaaEVqubQAAQHDz27BTt25dE1ji4+M9OhJrX5xWrVqZ2/ozJSVFNm/e7N5m1apVkp2dbfr2AAAA+LQZS6+Hs2fPHo9OyVu3bjV9bmrVqiUjRoyQiRMnSoMGDUz4GTdunBlh1aNHD7N9w4YNpUuXLvLAAw+Y4ekZGRkybNgwM1KroJFYAAAguPg07Hz77bdy8803u2+PHDnS/Ozfv7/MnTtXHn/8cXMtHr1ujtbgtGnTxgw1j4iIcD/m3XffNQGnQ4cOZhRWr169zLV5ABSPjMws2bdvnwQSHbVJ0zYQvHwadtq1a2eup1MQvary+PHjzVIQrQWaN29eMZUQQE6H0zIkMSlZ4oYPlrCwwBnfEFGuoixcuozAAwSpwPlrBcDn0tIzJSwkW55tX1ZiYyIlECQcSZdx8cdN7TBhBwhOhB0A561O5XCJjS4jgSPN1wUA4EN+OxoLAADAGwg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFiNsAMAAKxG2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFgt1NcFAIDilpGZJfv27ZNAEhUVJdHR0b4uBmAFwg4Aqx1Oy5DEpGSJGz5YwsIC509eRLmKsnDpMgIP4AWB85sPAEWQlp4pYSHZ8mz7shIbEymBIOFIuoyLPy4pKSmEHcALCDsAgkKdyuESG11GAkearwsAWIMOygAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNb8OO88884yUKFHCY4mNjXXfn56eLkOHDpXKlStLuXLlpFevXnLo0CGflhkAAPgXvw476sorr5TExET3sm7dOvd9jzzyiHz00UeyYMECWbNmjfz+++/Ss2dPn5YXAAD4F7+fCDQ0NDTfWX9PnDghs2fPlnnz5kn79u3Nujlz5kjDhg1l48aNcv311/ugtAAAwN/4fc3O7t27JSYmRi677DLp27evHDhwwKzfvHmzZGRkSMeOHd3bahNXrVq1ZMOGDYU+5+nTpyU1NdVjAQAAdvLrsNOyZUuZO3euLF++XGbOnCkJCQly4403SlpamiQlJUmpUqUkKirK4zHVqlUz9xVm8uTJUqFCBfdSs2bNYt4TAADgK37djNW1a1f3/5s0aWLCT+3ateX999+X0qVLF/l54+LiZOTIke7bWrND4AEAwE5+XbOTm9biXH755bJnzx7Tj+fMmTOSkpLisY2Oxsqvj09O4eHhEhkZ6bEAAAA7BVTYOXnypOzdu1eqV68uzZs3l7CwMImPj3ffv2vXLtOnp1WrVj4tJwAA8B9+3Yz12GOPSffu3U3TlQ4rf/rpp6VkyZLSp08f09dm0KBBpjmqUqVKpnbmoYceMkGHkVgAACAgws6vv/5qgs3Ro0flkksukTZt2phh5fp/NXXqVAkJCTEXE9QRVp07d5bXXnvN18UGAAB+xK/DznvvvVfo/RERETJjxgyzAAAABHyfHQAAgPNF2AEAAFYj7AAAAKsRdgAAgNUIOwAAwGp+PRoLAIJVRmaW7Nu3TwLtKvdnu4I94AuEHQDwM4fTMiQxKVnihg+WsLDA+TMdUa6iLFy6jMADvxM4v0UAECTS0jMlLCRbnm1fVmJjAmPuvoQj6TIu/riZr5CwA39D2AEAP1WncrjERpeRwJHm6wIA+aKDMgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVCDsAAMBqhB0AAGA1wg4AALAaYQcAAFgt1NcFAADYISMzS/bt2yeBJCoqSqKjo31dDBQzwg4A4IIdTsuQxKRkiRs+WMLCAuejJaJcRVm4dBmBx3KBc0YCAPxWWnqmhIVky7Pty0psTKQEgoQj6TIu/rikpKQQdixH2AEAeE2dyuESG11GAkearwuAi4AOygAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYAQAAViPsAAAAqxF2AACA1Qg7AADAaoQdAABgNaaLAAAELWZqDw6EHQBAUGKm9uAROEcXAAAvYqb24EHYAQAENWZqtx8dlAEAgNWo2QEAIIDQqfr8EXYAAAgQdKoumsB5pwAACHJ0qi4awg4AAAGGTtXnhw7KAADAaoQdAABgNcIOAACwGmEHAABYzZqwM2PGDKlTp45ERERIy5Yt5ZtvvvF1kQAAgB+wIuzMnz9fRo4cKU8//bRs2bJFmjZtKp07d5bk5GRfFw0AAPiYFWHn5ZdflgceeEDuu+8+adSokcyaNUvKlCkjb731lq+LBgAAfCzgr7Nz5swZ2bx5s8TFxbnXhYSESMeOHWXDhg35Pub06dNmcTlx4oT5mZqa6tWynTx5UrKysmXH73/IyfQsCQS7D/8p2Y4jOxNPSaYTGKcHZb54ArHclPnioMwXRyCW+Zdj6eazUD8Tvf0563o+x3EK39AJcL/99pvuobN+/XqP9aNGjXKuu+66fB/z9NNPm8ewsLCwsLCwSMAvBw8eLDQrBEYs9DKtBdI+Pi7Z2dly7NgxqVy5spQoUaLA9FizZk05ePCgREYGxiW6vYV9Z9/Z9+ARzPse7PufGoD7rjU6aWlpEhMTU+h2AR92qlSpIiVLlpRDhw55rNfbBc3BER4ebpbcM7KeCz0BAuUk8Db2nX0PNux7cO57sO9/ZIDte4UKFezvoFyqVClp3ry5xMfHe9TU6O1WrVr5tGwAAMD3Ar5mR2mTVP/+/eXaa6+V6667TqZNmyZ//PGHGZ0FAACCmxVh56677pLDhw/LU089JUlJSdKsWTNZvny5VKtWzWuvoc1eeh2f3M1fwYB9Z9+DDfsenPse7PsfbvG+l9Beyr4uBAAAQHEJ+D47AAAAhSHsAAAAqxF2AACA1Qg7AADAaoSdczBjxgypU6eORERESMuWLeWbb76RYLB27Vrp3r27uTKlXll68eLFEiwmT54sLVq0kPLly0vVqlWlR48esmvXLgkGM2fOlCZNmrgvLKbXq1q2bJkEo+eee86c+yNGjBDbPfPMM2Zfcy6xsbESLH777Te55557zJX0S5cuLY0bN5Zvv/1WbFenTp08x12XoUOHik0IO2cxf/58cx0fHY63ZcsWadq0qXTu3FmSk5PFdnqtIt1fDXvBZs2aNeaXfePGjbJixQrJyMiQTp06mffEdjVq1DAf8jrBrv6xb9++vdx+++2yY8cOCSabNm2S119/3QS/YHHllVdKYmKie1m3bp0Eg+PHj0vr1q0lLCzMBPudO3fKSy+9JBUrVpRgOM8Tcxxz/Xun/vrXv4pVvDkpp410MtGhQ4e6b2dlZTkxMTHO5MmTnWCip8qiRYucYJWcnGzegzVr1jjBqGLFis6bb77pBIu0tDSnQYMGzooVK5ybbrrJGT58uGM7nSC5adOmTjAaPXq006ZNG18Xwy8MHz7cqVevnpOdne3YhJqdQpw5c8Z8u+3YsaN7XUhIiLm9YcMGn5YNF9eJEyfMz0qVKkkwycrKkvfee8/UaAXT9Ctaq3fLLbd4/O4Hg927d5tm68suu0z69u0rBw4ckGCwdOlScwV+rc3QZuurr75a/vWvf0kwfub95z//kYEDBxY4KXagIuwU4siRI+aPfe4rMettvVIzgoPOtaZ9NrSa+6qrrpJgsG3bNilXrpy5kuqDDz4oixYtkkaNGkkw0HCnTdbabyuYaH/EuXPnmqvPa7+thIQEufHGG82M0rbbt2+f2ecGDRrIZ599JkOGDJGHH35Y3n77bQkmixcvlpSUFBkwYIDYxorpIoDi/pa/ffv2oOm/oK644grZunWrqdFauHChmXtO+zHZHngOHjwow4cPN/0WdEBCMOnatav7/9pPScNP7dq15f3335dBgwaJ7V9otGZn0qRJ5rbW7Ojv/KxZs8y5Hyxmz55tzgOt3bMNNTuFqFKlipQsWVIOHTrksV5vR0dH+6xcuHiGDRsmH3/8saxevdp03A0WpUqVkvr160vz5s1NDYd2VP/nP/8pttNmax18cM0110hoaKhZNORNnz7d/F9reoNFVFSUXH755bJnzx6xXfXq1fME+YYNGwZNM5765ZdfZOXKlXL//feLjQg7Z/mDr3/s4+PjPb4B6O1g6r8QjLRPtgYdbb5ZtWqV1K1bV4KZnvenT58W23Xo0ME04WmtlmvRb/zaf0X/r19+gsXJkydl7969JgjYTpuoc19a4ueffzY1W8Fizpw5pr+S9lWzEc1YZ6HDzrUaU//gXXfddTJt2jTTWfO+++6TYPhjl/Nbnbbh6x987aRbq1Ytsb3pat68ebJkyRJzrR1XH60KFSqYa3DYLC4uzlRl6zHW/hr6PnzxxRemL4Pt9Fjn7pdVtmxZc+0V2/trPfbYY+a6WvoB//vvv5vLbWi469Onj9jukUcekRtuuME0Y915553mWmpvvPGGWYLly8ycOXPMZ53WYFrJ18PBAsErr7zi1KpVyylVqpQZir5x40YnGKxevdoMt8699O/f37Fdfvuty5w5cxzbDRw40Kldu7Y53y+55BKnQ4cOzueff+4Eq2AZen7XXXc51atXN8f90ksvNbf37NnjBIuPPvrIueqqq5zw8HAnNjbWeeONN5xg8dlnn5m/b7t27XJsVUL/8XXgAgAAKC702QEAAFYj7AAAAKsRdgAAgNUIOwAAwGqEHQAAYDXCDgAAsBphBwAAWI2wAwAArEbYASyk1wodPHiwmdqjRIkSZpoPFO7o0aNmbqD9+/f7uijwgZ07d5rJfnU6INiHsIOAM2DAAPMBrktYWJiZpPPxxx+X9PR0j+1c22zcuNFjvU5oqXMd6X0655OLzm7dvn17ExDKlCkjDRo0MHPFnDlzJt9y9O7dW7p06eKxbvny5eZ5n3nmGY/1evtiziem5Zg7d66ZsT0xMbHY53XSfV68ePE5bZff8t5774mv/eMf/5Dbb79d6tSpI8FCzxGd3fxc6O+KzgYfHh4u9evXN48tjIbG/I51zt9HfY7c90dERIgv6Kzn119/vbz88ss+eX0UL8IOApKGDP0Q37dvn0ydOlVef/11M3FhbjVr1jQT3OWkM5mXK1cuz7c6fU6d8HXt2rVm5utXXnnFzHyflZWVbxluvvlm+eqrryQzM9O9bvXq1eY1c4Yo13rd/mJxzVatkxtGR0fnO7lfQSGuuOnx0GOXc+nRo0e+2+p7r5MUeqvsBT3u1KlTMnv2bBk0aFCRntd2Ogmwzoat57DWEo4YMULuv//+c5ocduXKlR7Hunnz5h73R0ZGetz/yy+/iK/oBM8zZ870+J2GJXw9ORdwvnQi0ttvv91jXc+ePZ2rr77aY52e3mPHjnUiIyOdU6dOudf/5S9/ccaNG2fu18lO1dSpU506deqcVzl00jx9jg0bNrjX6USxM2bMcCIiIpw///zTrNOfOrmgaxLRxx9/3GnQoIFTunRpp27duqaMZ86c8XjOH3/80eO1Xn75Zeeyyy5z3962bZvTpUsXp2zZsk7VqlWde+65xzl8+LD7/ck5ealO6uma0HLo0KFmUsvKlSs77dq1M+u/+OILp0WLFmYCyOjoaGf06NFORkaG+7X0cQ899JAzatQop2LFik61atWcp59+2n2/Pn9+r5cfvX/RokUF3q/vUYUKFZwlS5Y4DRs2dEqWLOkkJCSY5xw/frxz7733OuXLl3dPRrtw4UKnUaNGpuy6zYsvvujxfAU9LrcFCxaYSU9zat68ufPCCy+4b+s5Fxoa6qSlpZnbBw8eNPuze/duc/udd94xjylXrpx5j/r06eMcOnTI3JeVlWUm13zttdc8XmPLli1OiRIlnP3795vbx48fdwYNGuRUqVLFlPfmm292tm7d6t5e/6/HTV9D77/mmmucTZs2Ffh+vvTSS2ZyyzJlyjg1atRwhgwZ4i5/fhP95jyuOek5e+WVV3qs04lCO3fuXOBr63HT5/zuu+/Oerwv9Pdfz2k9T130/8OGDTPro6KizO+ITux58uRJZ8CAAeb9q1evnvPpp596PM/p06fN7+rKlSvPq0zwf9TsIOBt375d1q9fb2phctNvkdos8cEHH5jbBw4cMDU39957r8d2Wvuh3yr1vnN1+eWXS0xMjKm1UWlpabJlyxb561//al5zw4YNZr2WTZvOXDU75cuXN9X3Wpv0z3/+U/71r3+Z2inXc2rt0rvvvuvxWnr77rvvNv9PSUkxzW1XX321fPvtt6bJ6tChQ3LnnXea+/U5x48fb/of6D5t2rTJ/Txvv/22eZ+0RmrWrFny22+/Sbdu3aRFixby/fffm2+1WsMxceJEj9fXx5UtW1a+/vprmTJlinn+FStWmPtcz++qscn5ekWhtSzPP/+8vPnmm7Jjxw7Tj0a9+OKL0rRpU/nuu+9k3LhxsnnzZrPP2pyoNXHaVKjrczev5H5cfr788ss8NQ433XSTu4ZOc5puo00+69atczd7XnrppaZJR2VkZMiECRPM+6hNetqMo02uKiQkRPr06SPz5s3Lc1xbt24ttWvXNrf13ElOTpZly5aZ/dNmow4dOsixY8fM/X379jXHVd9jvf+JJ54wTbkF0dedPn26eR/1GK5atco0+Sqt9Zs2bZpHzcpjjz2W7/PoudyxY0ePdZ07d3af44W57bbbzDFs06aNLF26NM/9J0+eNPuvNaLajKhl9Qbd3ypVqsg333wjDz30kAwZMsS8v7rf+nvaqVMn83dAzzcX/d1o1qyZOdawjK/TFnC+9JudfuPXWg39FqancUhIiPmWn18twrRp08w3ZPXss886d9xxh/kGnbNmJzMz03zj03Vau9GjRw/nlVdecU6cOFFoWfr27et06tTJ/P+TTz4xtQxq8ODBzlNPPWX+r7VIWoNTEK090BoBF61l0m+dLrlreyZMmOB+TRdXLYNu63qO3DUs+m03d+3XmDFjnCuuuMLJzs52r9OaKf3mq7URrse1adPG43FaE6Q1QLnf67PR7bTWS49dzuWXX35xf9PXbXLWZijdFz0mOd19992mli4nrX1yHYOCHpcfrSkYOHCgx7qlS5eaWgc9N7Q8el5oTYFrv++//35ThoJojYvui6smRWs4tBbHta+u2p6ZM2ea219++aWphUxPT/d4Hj0XXn/9dfN/rc2ZO3euU1Rag6W1eudbs6I1kZMmTfJYp+e77l/OWtOctKZRa5Y2btzofPPNN+Z90/3XWjuX9evXO2+//bZ5b7SG8dZbbzXvgZ7PF1qzk/Oc1WOo55nW8LkkJibmqZlV+vdB/xbALtTsICC5+g5oTYN2Ita29l69euW77T333GO+gWr/Hv3WP3DgwDzblCxZ0tRM/Prrr6bmQr+xT5o0Sa688krzjbcg7dq1M7Uk+q1eawH0du5aAf2Zs7/O/Pnzzbd5rU3SvkNjx441NU4uWlOhtQKujpz67V+/4cfGxprbWnOgtUn6WNfiuk/76hQmd+3Fjz/+KK1atTIdQ120bPptW98LlyZNmng8TvsDaQ1EUWgtlh67nIvWkOX8dp379ZTWeOUuu5Y1J729e/duj35WuR+Xnz///DNPx9gbb7zR1NZpjZDW4ugx1ePrOq66znW8lda0dO/e3XRE19o73V65jq3WGDRs2NBdu6OP1/dQaxtcx1Xfd+08n/PYan8Z13EdOXKk6SujtSzPPffcWY+39pfRmiE9n7VMWpOho85y1mYUF61V0fK2bNnS1BxqefV38YUXXnBvo+dev379zHuj79eHH34ol1xyiemDd6FynkP6+63va+PGjd3rqlWrZn7mPo9Lly59Ud4fXFyEHQQkbVLR5gNtnnjrrbdM6NHml/zoH7lbb73VdD7VEVtdu3Yt8Hn1Q0E/EF599VVTna7ba3NPQTTE6FBVbVbQAOL6gNOfWiZtftCf2uykNHRpU4Q2HelIKf0gffLJJz06zmoI0u1dH4r6Ux/joh+I+qGaOzDoh3zbtm3P+r4VRe6mEg1H+XUcPhe6f3rsci45O1Drh03O8HWhZT+Xx+kH8/Hjxz3WaZOVnl8ablzBRt9fPWY///yzeb9dx1vPAW3W0SYhDad6PmhHeJXz2OpxzHlctVO8np+u46ohMvdx3bVrl4waNcpso011el5qZ2FtktIRRK7XyU0Ds573+qGvzbgaxmbMmJGnTOd6zLSpNCe9rfurx+tcafDZs2dPoeeZNs8Wtk1+8htEkN85m3Od6xzLfR7r76wGLtiFsIOAp/0SxowZY2pI9Bt6frQ2Rz+09Fukfss7FxUrVjQfPoVdd6NevXqmr4H2RdAPJteHn4YmXV566SXzweKq2dH+O9o/QQOO1jjo8Pb8Rp/oh6LWALlqpLS2x0VrefQDT/sF5Q4N5xsItKZBX+N/LUz/ozVVWgugfUPOlX6IFDRqrbho2bWsOelt7fd0rsfYRT9gtQ9Vbno8NcRqXy4NO3pZAn1dHaau54a+lvrpp59MjYnWXmiNkNa05Vfzpf2utI+ZBo+FCxd6hFg9rklJSSb45T6uGsZc9DUfeeQR+fzzz6Vnz555Rhu66GvoB7megzqkWh/3+++/e2xT2GjDnLQGJj4+3mOd9tnS9edDf0f0fSuIlkX7XxW2jcodvPR3xFv0+Oj5ALsQdmAFbQrQDzjXN9fc9Bv04cOHTcfa/Gi1uXZg1A8QbRrQMDF69GjzU2tRCqNB5rXXXjMfSq6qcdcHpQ5fd3VkVhputFlDryujr6OdR/P7Zq4fYtqEomXS58/ZzDN06FDz7VM7vGoNgj6PDgHWprzzDRx///vf5eDBg6YDp35gL1myxAzh1+YHDZHnSoOXfhjqh3XuGpLctIO1bpdzKcqF3B599FHzmtopWGtatEOq1sgV1Mm2MForo8c6d9k14Oh7qwHE1VSo67T2xhVslTZdaXDQ460fvBp+tVz5vU/aQVZrGfVYaeddF22a0vCgw/D1PNSaGQ3HGoy1I7oG+WHDhpnQrgFZg50efw1f+dHzUZtXXWX697//naeWUsujNUr6Ph45cqTA5psHH3zQPId2btbzRM/3999/34QuF33vtcnMRY/Hf//7X7O9LtosrLWweq656O+j7qs+t3Ya1mYu3TdtqiuM1pZqx359nHZk12Okzc3a5Hch9D3XTvu5O2PDAr7uNAScr/w6KKrJkyeb4cM6vPRsnWZzd1DWIcA6fFs7EmunZ+3E2bZtW9NJ9WxcnWoffPBBj/XakVTX/+1vf8vTiVafXzsB6/Bd7UycXyfRO++80zz+rbfeynPfzz//bDpS6rBaHcIeGxvrjBgxwt3RuKAOytqRM7dzGXqe+3H6/uccxq3vU/369c3Q7LMNPc9v0WPnei/zey/0OXWfcnMNPQ8LC3Nq1arlMVS8sMflRy8bMGvWLI91R48eNZ1q9Ti56DmlZc697bx588zlC/T8adWqlXlP8ht6rcPPdX2/fv3ylCE1NdUM84+JiTH7VLNmTdMJ/sCBA2ZYdO/evc06PVa6jQ6vdl3iID96yYLq1aubc0SHievweH1tPf9d9LzV87GwoedKf1eaNWtmXlsvg+C6lIKLPjbnsdfzXy8foMPetdOxvr/aQTonPWf1uOlz6nD9bt26md/Fwuh51759e7M/+jh9Xt0v7bytQ+sLOmfzOxdy/43QTtiFDadH4Cqh//g6cAGAr33yySemb4w2Y5xPrRYuLh3Or7WD53LF7vOhzc1a86p9qXJ3fEfgy3tZVQAIQtrpVzsdazOG9sNCcNHmZe37R9CxE2EHAP4/nQYBwcnVGRx2ohkLAABYjYZpAABgNcIOAACwGmEHAABYjbADAACsRtgBAABWI+wAAACrEXYAAIDVCDsAAEBs9n9I+oZGlr4KuAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(rms_wavefront_error, color=\"C1\", edgecolor=\"k\", alpha=0.8)\n",
    "plt.xlabel(\"RMS Wavefront Error (waves at 0.55 µm)\")\n",
    "plt.ylabel(\"Occurrences\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's dive deeper into the relationships by looking at the correlation matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns  # pip install seaborn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame(\n",
    "    dict(\n",
    "        radius_of_curvature=radius_of_curvature,\n",
    "        refractive_index=refractive_index,\n",
    "        rms_spot_radius=rms_spot_radius,\n",
    "        rms_wavefront_error=rms_wavefront_error,\n",
    "    ),\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAd4AAAGUCAYAAACbXUPJAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAATSZJREFUeJzt3QeYFNXyNvBagqCg5BwkSs5JkoAgUQQBkaSASxYRQRD+lxwkBwmCCIgoqEgSQXJOAiKoKDl6kaQkAYk73/PWteebGTYNLL3bp9/ffeZxt2d2p2eWO9VVp845IR6PxyNERERki3j2PA0REREBAy8REZGNGHiJiIhsxMBLRERkIwZeIiIiGzHwEhER2YiBl4iIyEYMvERERDZKYOeTERERHZ03I6jH52wcKiZh4CXH/5/SafAhcubPS2K6DKlTSJX+Zv8tYf3AUBm1eKOYrkf9SrF9CsZg4CUiIluFhISImzHwEhGRvULc3V7k7ldPRERkM2a8RERkr3gsNRMREdkmxOWlZgZeIiKyV4i7M153X3YQERHZjBkvERHZKiRefHEzBl4iIrJXPJaaiYiIyCbMeImIyFYh7GomIiKyUYi7S80MvEREZK8Qd2e87n71RERENmPGS0REtgphqZmIiMhG8dxdbGXgJSIiW4Uw4yUiIrJRCAMvERGRfULcXWp296snIiLjTJ48WbJlyyaJEyeWMmXKyM6dOyN87McffywVK1aUFClS6K1atWqRPj4mMPASEZGtQuKFBHULxldffSXdunWT/v37y48//ihFihSRGjVqyPnz58N9/IYNG6Rp06ayfv162b59u2TJkkWqV68up0+fFlcE3lmzZkny5Mm93w8YMECKFi0qccXWrVulUKFCkjBhQqlfv35snw4RkXNLzSFB3IIwduxYadu2rbRu3Vry588vU6dOlSeeeEJmzpwZ7uPnzJkjnTp10liTN29emT59uoSFhcnatWvFFYE30LvvvvtIX3ywcBWFP87x48f1IsGpcIWHrsLLly/H9qkQkVubq0Kif7t165ZcvXrV74ZjgW7fvi27d+/WcrElXrx4+j2y2ei4ceOG3LlzR1KmTCmOCrx48TEhadKkkipVqhj5XTHh6NGj8vzzz0vmzJn9MnO7xNT7GpPwD5SI6FEaNmyYJEuWzO+GY4H+/PNPuXfvnqRLl87vOL4/e/ZstJ7rvffek4wZM/oF7zgZeCtXriydO3eWrl27SurUqbWejnQfZdkkSZJozRyp/LVr1/x+Dllj1qxZtQzw8ssvy19//eV3f2CpGc+D5/CFkm+rVq2833/44YeSO3duHVTHm92oUaNovQZcPXXp0kXSpk2rP1uhQgXZtWuX3nfixAnNEHF+b7zxhn4dnYz3119/lRdffFGeeuopefLJJ3UAH8E7uq8FzQGDBw+W119/XX9Hu3btpFy5cvoPw9eFCxe0/L1p0yb9/rPPPpOSJUvqc6ZPn16aNWvmHd/Aa6lSpYp+jUYCvBbrOfF848eP9/vdeP/xd7Dg8VOmTJGXXnpJ/7ZDhw7V4998840UL15c37scOXLIwIED5e7du9F674nIfbsThQRx6927t1y5csXvhmMxbfjw4fLll1/KokWL9LMszme8n376qTz22GM6DoqaOtL7CRMmaPDBfevWrZOePXt6H79jxw4JDQ3VgL13714NBkOGDHmoc/jhhx80eA4aNEgOHjwoK1askOeeey5aP4tzW7BggZ4rBuRz5cqlFxAXL17UC4czZ85o8ENgwtevvvpqpL8PA/N47kSJEulrR/kDQTvYYDR69GhtDtizZ4/07dtXmjdvrv8wPB6PXzMBrtAQ2K0sFAH7p59+ksWLF2uwtYIrXgteJ+A9wmv54IMPgjonBGJcKP3yyy/6mjZv3qwXB2+//bb89ttv8tFHH+mFiRWUiYh8hcSPF9QNn6P4/PW94VggJH7x48eXc+fO+R3H90hCovqsReBdtWqVFC5cWBwxjxdZ5siRI73f58mTx/s1MikE1Q4dOmhGCviwr1mzpjcYP/PMM7Jt2zYNlg/q1KlTmoUhy0S29/TTT0uxYsWi/Lnr169rFodgUatWLW+L+erVq2XGjBnSo0cP/aMh20OJI6o/oNXOjsciSCIbtV5jsFDa7t69u/f7xo0ba6a8ZcsWb6CdO3euduVZq8EgGFqQfeICqFSpUlpxQPneGrtAdv8gJXNk0GhcsOD5evXqJS1btvQ+JwI//rboLIysyhA4ThPe/5mIiKIDyV+JEiW0N8hqgLUapZDkRQSxC4nCypUrtVr4qMVYxosX62vNmjVStWpVyZQpkwbB1157TUu1GLiG/fv36/wqX2XLln2oc3jhhRc02OKDH8+HbjXr+SKD8i+yxPLly3uPIViWLl1az/NBIItHYLSC7oMK/EeQJk0abXXHawM0eqFpAJmwBdl13bp1tYyP975SpUreC5OYEHhOyKxRZUBQt27oKkQ2Hdn7H91xGyIyTMij62pGEywSJ1Qv8fndsWNHTa6sZAHVOd8y9YgRI7SaiK5nJIkYC8YtcGg0TgZeZJoWlDaRdSJdR1kTgQAZ4MM2CKF87VtiDWzuQZBBmfiLL76QDBkySL9+/bRMGxvdu48//vhDvZbw3lcLguz8+fP18ch2MZaOG+AfGErkKMUgOGOcGuMV0XnvH/Sc8A8UY7q42LBuKEMfPnw40nESu8ZtiMjZXc3BwDAgysb4/EePCj6PUEm1Gq6QgCApsKDaic9G9AMhblg3/A5HLRmJQIv0fsyYMfphDvPmzfN7TL58+XSc19f3338f6e9Ftuf7hqF7bd++fd5mIUiQIIF2o+GGMidKqRhjbdCgQYS/N2fOnN7xaWTMVsBB0ApsgIouXHTgigu/J7ysNzqvJSL16tXTRiv8Y0LgxRWc5cCBA1pZwFgFxnOtsW9feK3Wc0Z2TmjZR0YdFTRVYbwY4+LBQFmZpWUi9wl5xEtGoqwcUWkZ0yl9IVG02yN59fgARsCZOHGiHDt2TLts0XDlC01QCBy4qkBmNGnSpCjHdzHeuWzZMr0hwKCE4JvNLl26VMczcYVz8uRJmT17tl4A+I43hwcZHH4XxnJxDmgQQqkUZVI0gD0I/NERuJo0aaKBD68R7wMCVHReS1Tni/ELlEdQSsH4rgXlZQRW671fsmSJjrf6wsUFxoPxfqEj2iqp4JxwjmiWQsaKMVs0KkQFV5Z4r5H1opkO54Sx7T59+gT5rhERme+RBF6UdzGdCLXzggULaskzcOzu2Wef1To8mqzweHSSRfVBjSYeBANkeBi3xFiub4aI7HbhwoUaQJBRI9ij7FygQIEozxkZYsOGDXVsGBnckSNHdKAdU24eBOYfI9NGUMO5Ygwcr9fKfqN6LVFBuRljqxhHRrD1zVrRJPb111/rqi14XYElE4y7I0iiIQrlF+vKEGVenAuGCerUqaPBHdWAqKC0jSCOvyGauPC3HTdunLd6QETkJ15IcDfDhHgCB/WIYtnReTPEZDkbh8qZPy+J6TKkTiFV+pv9t4T1A0Nl1OKNYroe9f/XpBkTzuzaEtTjM5SqICbhtoBERGSrEJfvxxun12qOKehi853qEnh7kGk2mJMc0e/DfUREZP90IidwRcaLVZ3QcBXZ/cHCvFVs4hAeTOUhIqLwhbg843VF4MUUo2CnukQFqz7hRkREFAxXBF4iIopD4jHjJSIisk+IeeO2wWDgJSIiW4W4fIzX3ZcdRERENmPGS0RE9gpxd87HwEtERPaK5+5SMwMvEREZtTtRXMfAS0RE9gphxktERGSbEAZeIiIiG4W4u9Ts7ldPRERkM2a8RERkr3gsNRMREdkmxOWlZgZeIiKyV4i7M153X3YQERHZjBkvERHZK8TdOR8DLxER2SrE5aVmBl4iIrJXPGa8REREtglxecYb4vF4PLF9EkRE5B6Xfj8R1ONTZMkmJmHGS3HOmT8vxfYpPFIZUqeQo/NmiOlyNg6V44s/F9Nlr99Czu7eLqZLX6JsbJ+CMRh4iYjIXiEc4yUiIrJNCJeMJCIislGIuzNed796IiIimzHjJSIiW4XEiy9uxsBLRET2CnH3GC9LzURERDZixktERLYKYVczERGRjULcXWxl4CUiInuFuDvjdfdlBxERGWfy5MmSLVs2SZw4sZQpU0Z27twZ6eO//vpryZs3rz6+UKFC8t133z3S82PgJSIiW4WExAvqFoyvvvpKunXrJv3795cff/xRihQpIjVq1JDz58+H+/ht27ZJ06ZNJTQ0VPbs2SP169fX2759++RRYeAlIiJ7xQsJ7haEsWPHStu2baV169aSP39+mTp1qjzxxBMyc+bMcB//wQcfSM2aNaVHjx6SL18+GTx4sBQvXlwmTZokjwoDLxER2SskXlC3W7duydWrV/1uOBbo9u3bsnv3bqlWrZr3WLx48fT77dvD30EKx30fD8iQI3p8TGDgJSIiW3lCQoK6DRs2TJIlS+Z3w7FAf/75p9y7d0/SpUvndxzfnz17NtxzwfFgHh8T2NVMRES2uhcW3ON79+6t47a+EiVKJE7FwEtERHFaokSJohVoU6dOLfHjx5dz5875Hcf36dOnD/dncDyYx8cElpqJiMhWniD/F12PPfaYlChRQtauXes9FhYWpt+XLVs23J/Bcd/Hw+rVqyN8fExgxktERLbyRD+WBg0l6ZYtW0rJkiWldOnSMn78eLl+/bp2OcPrr78umTJl8o4Rv/3221KpUiUZM2aM1KlTR7788kv54YcfZNq0aY/sHBl4iYjIVmGPMPK++uqrcuHCBenXr582SBUtWlRWrFjhbaA6deqUdjpbypUrJ3PnzpU+ffrI//3f/0nu3Lll8eLFUrBgwUd2jiEez6O89iAK3pk/L4nJMqROIUfnzRDT5WwcKscXfy6my16/hZzd/eimnsQV6UvEXOn1wqUrQT0+TYpkYhLjx3hxXdGuXTtJmTKlhISEyN69e219/g0bNujzXr582Zbna9Wqla668rBwzrjqIyJ6FJ/LniBupjG+1IwSw6xZszQA5siRQ7veHpXKlStrWQNjCr5ljDNnzui8MztgFRYT/6ESkTnCXP4Z5ejAi1VK0MUWmaNHj0qGDBk0AD7M73lQ+L2Psi09kF0BnojoQXncHXedVWpGRtm5c2fp2rWrZq5Y1gsLWdeqVUuSJk2qg+evvfaarl5ilV3feustHUxH6RS7VUT0e6w1PrEzRZIkSSRLlizSqVMnuXbtmt85bN26VX8ea3+mSJFCf/bSpUv6XBs3btSME8+F24kTJ/xKzVjm7PHHH5fly5f7/c5FixbJk08+KTdu3NDvf//9d2ncuLEkT55cS+T16tXT3/UgpWaca5cuXaRnz576u3ARMGDAAL+fOXz4sDz33HO6MwfWNkUrfaDIzunAgQP6fqBBwTJv3jx9rb/99lu0zpuI3MPj8lKzowIvfPrpp5pFIgAOHz5cnn/+eSlWrJi2f6OsjInPCBCAIDho0CDJnDmzlnt37doV7u/BItqATrcJEybIr7/+qvevW7dOA5YF48NVq1bV4IR1PLds2SJ169bVJcrwXJj3hcW58Vy4IXj7euqpp+TFF1/0C1AwZ84cDZYIXnfu3NFgjkC8efNmPT9cVGARb2TmD/qe4WJix44dMnLkSH1PrOCKOW4NGjTQ9wL347147733/H4+qnPCdlqjR4/WCxVc5Pz3v/+VDh06yIgRI/S9IiIKLDWHBXEzjeNKzWj1RvCAIUOGaNB9//33vfdjBwoEvEOHDskzzzyjwQIrmQSWe31/jwUZsAXZMX4/AsiHH36ox/B4zA2zvocCBQp4v0bwQvCMrLTcvHlzzcqR3eKxyIKXLVumWa+1pRWC4fTp0zVThk8++UQzTWTP1atXD/o9K1y4sG6RZb1u7LqBCeMvvPCCrFmzRjPWlStXSsaMGfUxeD9RRbBE55wQdLGHZYsWLfR9KFWqlFYbIoNFzgMXOnfyMnBEFD0e82Kp2RkvViWx/PTTT7J+/XrNvqwbsi9rbDe6v8eCIISMFpOrEbARIP/66y9vCdjKeB9G7dq1JWHChLJkyRL9fsGCBZoJW7tj4DUdOXJEn996TSjt3rx5M8rXFFng9YUxb2tvyv379+uFihV0IXDFluieEy56fv75Z90DEw1tVpCOSHQXPiciMonjMl6UTC0Yf0WpFyXNQAgu0f09gPFKlIE7duwoQ4cO1cCCUjI2R0Y5FdkpxiwfFrLBRo0aabm5SZMm+l9M+E6QIIH3NeGiAOXnQGnSpHmg50Sg94WAiAw2uqJ7TgjQWCEGJXuU2qP6G0S08PnFv/93oUNEZvK4POV1XOD1hc2KkTGiLGwFrgeFPRwRjLBsmLWqCRqEAjNHlGgHDhwYYVDFeG9UUG5GmRdjyRhHRknb9zWhtJs2bVrNhB81bPyMxinfQPn999/7PSY653Tx4kVt7PrPf/6jvwuvEZlvZBcrES58zsBLZLQwlwdex5Wafb355pv6gd+0aVNtnELZE2OVWJMzOgHQV65cubSJaOLEiXLs2DH57LPPvE1XvhkangfjmSipYmx0ypQp3i5qXACgQQnZM45FlFWigxjjwAhO2bNnlzJlynjvwzF0WqNrGI1Mx48f13FUdCajaSmmocSNsXCsbYqMFc+J4OkrOueEsXCUrLHsGrrD8f6/++67MX6+RERO5+jAi3FJdNjiQx4NPpgKhAYpNP34rsUZHUWKFNGAgbI11uhEWTVwvBEBatWqVRqgsPg2xkK/+eYbb7aNQINGLnTyogSLDt/woNSLiwX8HgQ1Xyhpb9q0SbJmzardxshIUe7GeOqjyIDxPqGx659//tHX1KZNGy21B3NOs2fP1sYqXKzgvUAZ//PPP5ePP/74vqlTRESeIG+m4VrNFOdwrWYzcK1ms8TkWs1HTvvvfxuVXJn+t8GBKRyd8RIRETmNo5ur3AhTeSKCsm7FihVtPR8iomB5XF5nZeB1mMh2V8L8YyKiuC7M5ZGXgddh0H1NRORkHpcHXo7xEhER2YgZLxER2crj7oSXgZeIiOwV5vLIy8BLRES28jDwEhER2SfM3XGXgZeIiOzlMXIhyOhj4CUiIlt5XF5q5nQiIiIiGzHjJSIiW4W5O+Fl4CUiInt5XF5qZuAlIiJbeVweeDnGS0REZCNmvEREZKswl2e8DLxERGQrj7vjLgMvERHZK8zlkZeBl4iIbOVh4CUiIrKPx91xl13NREREdmLGS0REtgpzecob4nF7sZ2IiGy1ad+RoB7/XMFcYhJmvBTnVOk/Q0y2fmCoHF/8uZgue/0WcnSe2X9LyNk41DWvM6Z4XJ7ucYyXiIjIRgy8RERk+xhvWBC3R+XixYvSvHlzeeqppyR58uQSGhoq165di/Txb731luTJk0cef/xxyZo1q3Tp0kWuXLkS1POy1ExERLa6F0f2BUTQPXPmjKxevVru3LkjrVu3lnbt2sncuXPDffwff/yht9GjR0v+/Pnl5MmT0qFDBz02f/78aD8vAy8REbnO/v37ZcWKFbJr1y4pWbKkHps4caLUrl1bA2vGjBnv+5mCBQvKggULvN/nzJlThg4dKi1atJC7d+9KggTRC6ksNRMRka08Hk9Qt1u3bsnVq1f9bjj2MLZv367lZSvoQrVq1SRevHiyY8eOaP8elJlRqo5u0AUGXiIiitOBd9iwYZIsWTK/G449jLNnz0ratGn9jiF4pkyZUu+Ljj///FMGDx6s5elgMPASEZGtwjzB3Xr37q2Zpe8Nx8LTq1cvCQkJifR24MCBh34NyLrr1KmjY70DBgwI6mc5xktERHFaokSJ9BYd3bt3l1atWkX6mBw5ckj69Onl/PnzfscxTovOZdwXmb///ltq1qwpTz75pCxatEgSJkwowWDgJSIiW3ke4RShNGnS6C0qZcuWlcuXL8vu3bulRIkSemzdunUSFhYmZcqUiTTTrVGjhl4ILFmyRBInThz0ObLUTERErpvHmy9fPs1a27ZtKzt37pStW7dK586dpUmTJt6O5tOnT0vevHn1fivoVq9eXa5fvy4zZszQ7zEejNu9e/ei/dzMeImIyJXmzJmjwbZq1arazdywYUOZMGGC937M7T148KDcuHFDv//xxx+9Hc+5cvmvH338+HHJli1btJ6XgZeIiFy5O1HKlCkjXCwDEEh9y+KVK1eOkTI5Ay8REbky8MYWjvESERHZiBkvERHZyuPuhJeBl4iIzJlO5AQMvEREZKswlwdejvESERHZiBkvERHZyuPyjJeBl4iIbBXm7rjLUjP9z4YNG3TXDqxdCrNmzdK9KomIYntbQNMw8MYRJ06c0MC3d+9eiQteffVVOXToUGyfBhEZyOPywOvYUvPt27flsccei+3TMPY9efzxx/VGRBTTwsTdHJPxYo1MLGbdtWtXSZ06tW7JhAxx5cqVUqxYMQ0Szz//vO6vuHz5ct154qmnnpJmzZp5F7iG+fPnS6FChfTxqVKlkmrVqulOE9EpxZYuXVqSJEmiJdjy5cvLyZMn9T5sgly0aFH56KOPJEuWLPLEE09I48aNdbNmC7aaGjRokGTOnFnPHY9fsWKF9/7s2bPrf/Fa8LrweqOCPSfr168vQ4cO1d008uTJo8c/++wzKVmypO4ViX0l8R4E7jv53XffyTPPPKPvQ5UqVTTj9hVYaraeyxf+Fr7n+aDvLRG5i8flGa9jAi98+umnmtFh+6apU6d6g96kSZNk27Zt8vvvv2vAGz9+vC58vWzZMlm1apVMnDhRH3vmzBlp2rSpvPHGG7J//34Npg0aNIjyD4vNkRF0KlWqJD///LNs375d2rVrpwHScuTIEZk3b558++23GlD37NkjnTp18t7/wQcfyJgxY2T06NH6O7Cf40svvSSHDx/W+61tp9asWaPnuXDhwmi9J2vXrtXdM1avXi1Lly717qgxePBg+emnn2Tx4sUaVH03hsb7hNddt25dLW23adNGevXqJQ/jQd9bIiK3cVSpOXfu3DJy5EjvBz0MGTJEs08IDQ2V3r17y9GjRyVHjhx6rFGjRrJ+/Xp577339GcQRBEQnn76ab0fGVpUsOcistcXX3xRcubMqceQUfu6efOmzJ49WzJlyqTfI9jXqVNHgy2yTgRcnAP2eoQRI0boeeEiYfLkyd6Nm5Ep4vHRhQx8+vTpfiVmBD8L3gdsc1WqVCm5du2aJE2aVKZMmaKvA+cGyJR/+eUXPacH9aDvLRG5j8flF+SOynhLlChx37HChQt7v06XLp2Wea2gax2zyqxFihTRfRcREF555RX5+OOP5dKlS9HaOgoZI7JUZInIXq3Ab8maNas36ELZsmW1vIxsFIH7jz/+8F4gWPA9ssOHgdcSOK67e/duPU+cE8rNyNTh1KlT+l88Z5kyZfx+Buf7MB7kvb1165a+N743HCMi86cThQVxM42jAi+yu0AJEyb0fo3Sr+/31jEEQIgfP76WZDEGnD9/fs1Kke1hA+OofPLJJ1piLleunHz11Vc6Pvr9999LXHtPMKaKCwSMb2OT5127dsmiRYu8zVcPCptEB16loqRteZD3dtiwYZIsWTK/G44REZnMUYE3JiAQI9McOHCgjsMiW7QCU1TQ+IRSNsaTCxYs6LeBMrJJZLUWBGUEKwQfBEE0P2Fs2he+R5ACK2u9d+/eQ72+AwcOyF9//SXDhw+XihUrSt68ee9rrEKZ3BpT9j3fyKAUHpjlB059Cva9xXuJEr7vDceIyGwelzdXOWqM92Ht2LFDm5GqV68uadOm1e8vXLhw33htIGRt06ZN02YoBFCUj9EU9frrr3sfkzhxYmnZsqWO5aJk2qVLF230ssZre/ToIf3799exVXQ0I4NG4EJWCjgfdAOjMQudz/h9yACDhfIyAh4yzg4dOsi+ffu00coXjmN8F+eExiqUptHFHBl0jI8aNUrHsVGW/vzzz/V342LkQd9bdHfjRkTucs/E+nEQXJXxIvPctGmT1K5dW0vFffr00QBUq1atSH8O48bIJBs2bKg/h47mN998U9q3b+99TK5cubSxCL8bwQdjzx9++KH3fgTibt26Sffu3XUcFAF2yZIl2jAGCRIk0CYoTElCcK9Xr94DvUZkpgiiX3/9tWbTyHxxMRAYnBcsWKAdzxibRYf4+++/H+nvRfm6b9++0rNnT23U+vvvv/0uPB70vSUicpsQj4l5vM0wpQlBLK6sOuV0VfrPEJOtHxgqxxd/LqbLXr+FHJ1n9t8ScjYOdc3rjCkzVvsPdUUl9IXSYhJXlZqJiCj2hbk832Pg/Rfmt0YEnbpoVLJbXDwnIqKH5XF33GXgtURWJvadnxtRqRm3uHROREQUNzHw+jRHxTVx8ZyIiB5WmMtTXgZeIiKyVRgDLxERkX08Lg+8rprHS0REFNuY8RIRka3uuXzlKgZeIiKylUfcHXhZaiYiIrIRM14iIrKVx90JLwMvERHZK8zlkZelZiIiIhsx4yUiIlt5XJ7xMvASEZGtwhh4iYiI7ONxd9xl4CUiInt5XB55GXiJiMhWYQy8RERE9vG4O+5yOhEREbnTxYsXpXnz5vLUU09J8uTJJTQ0VK5duxbtcnmtWrUkJCREFi9eHNTzMvASEZGtPB5PULdHBUH3119/ldWrV8vSpUtl06ZN0q5du2j97Pjx4zXoPgiWmomIKE6P8d66dUtvvhIlSqS3B7V//35ZsWKF7Nq1S0qWLKnHJk6cKLVr15bRo0dLxowZI/zZvXv3ypgxY+SHH36QDBkyBP3czHiJiMhWniBvw4YNk2TJkvndcOxhbN++XcvLVtCFatWqSbx48WTHjh0R/tyNGzekWbNmMnnyZEmfPv0DPTczXiIiitN69+4t3bp18zv2MNkunD17VtKmTet3LEGCBJIyZUq9LyLvvPOOlCtXTurVq/fAz83AS3HO+oGhYrrs9VuIG+RsbP7f0k2vM7ZKzYmCKCv36tVLRowYEWWZ+UEsWbJE1q1bJ3v27JGHwcBLcc6oxRvFZD3qV5Kzu7eL6dKXKCtH580QNwRdt7zOmBIW9ugaprp37y6tWrWK9DE5cuTQMvH58+f9jt+9e1c7nSMqISPoHj16VEvUvho2bCgVK1aUDRs2ROscGXiJiMgYadKk0VtUypYtK5cvX5bdu3dLiRIlvIE1LCxMypQpE2E23aZNG79jhQoVknHjxkndunWjfY4MvERE5LqVq/Llyyc1a9aUtm3bytSpU+XOnTvSuXNnadKkibej+fTp01K1alWZPXu2lC5dWjPh8LLhrFmzSvbs2aP93OxqJiIiV87jnTNnjuTNm1eDK6YRVahQQaZNm+a9H8H44MGD2skck5jxEhGRK6VMmVLmzp0b4f3ZsmWLMvA/yIUBAy8REdnKE/uV5ljFwEtERK4b441NDLxERGQrj8sDL5uriIiIbMSMl4iIbHXP5RkvAy8REdnK4/LAy1IzERGRjZjxEhGRrcLcnfAy8BIRkb08Li81M/ASEZGtwlweeDnGS0REZCNmvEREZCuPyzNeBl4iIrKVx91xl4GXiIjsFebyyMvAS0REtvIw8BIREdnH4+64y67mh3Xjxg1p2LChPPXUUxISEiKXL1+O7VMiIqI4jIH3IX366aeyefNm2bZtm5w5c0aSJUv2yJ7rxIkTGtz37t37yJ6DiMiOMd6wIG6meahS8+3bt+Wxxx4TNzt69Kjky5dPChYsaMz7dOfOHUmYMGGMvAanvXYievQ8Yl4wfWQZb+XKlaVz587StWtXSZ06tSRKlEgzsJUrV0qxYsXk8ccfl+eff17Onz8vy5cv14CEEmyzZs20JGuZP3++FCpUSB+fKlUqqVatmly/fj3S5963b5/EixdPLly4oN9fvHhRv2/SpIn3MUOGDJEKFSro1/fu3ZPQ0FDJnj27Pk+ePHnkgw8+8D521apVkjhx4vtKw2+//ba+BsuWLVukYsWK+juyZMkiXbp08Z4r3o8xY8bIpk2b9H3A95AtWzYZPHiwvP766/r627Vrp8cXLFggBQoU0PcNj8HP+sKx999/X9544w158sknJWvWrDJt2jTv/XgtgPfa9/miMn36dP1b4PXmzZtXPvzww/uy6K+++koqVaqkj5kzZ460atVK6tevL0OHDpWMGTPq+we//PKLvj/W3w6v7dq1a97fF9HPERH5rtUcFsRN3F5qRmkVGczWrVtl6tSpemzAgAEyadIkLbf+/vvv0rhxYxk/frzMnTtXli1bpkFu4sSJ+liUY5s2barBZf/+/bJhwwZp0KBBlF1uCFj4oN+4caN+j/Ku7/eAr61gFBYWJpkzZ5avv/5afvvtN+nXr5/83//9n8ybN0/vr1q1qiRPnlyDoQXBGgGoefPm3my2Zs2aOob7888/630IxLj4gIULF0rbtm2lbNmy+rrwvWX06NFSpEgR2bNnj/Tt21d2796t7wsuFBC88J7h+KxZs/xeJ4JxyZIl9ec6deokHTt2lIMHD+p9O3fu1P+uWbPmvueLCIIoXjsCId5vBHY8L/6Ovnr16qUXHXhMjRo19NjatWv1uVevXi1Lly7VCw7clyJFCtm1a5e+tzgX6/2wBP4cERE9RKk5d+7cMnLkSP0aH/5Wplm+fHn9Gllm7969NWjlyJFDjzVq1EjWr18v7733nv7M3bt3Ndg+/fTTej+y36ggK3vuuec0UOP34b+tW7fWbO7AgQOSM2dODfw9e/bUx6NUOnDgQL9scfv27Rp4EQDjx4+vQRAXBzhnK2AgA0aghWHDhmkQRoZvvfYJEyZoZjhlyhRJmTKlPPHEE3ohkj59er/zRVbYvXt37/f4PQj2CHrwzDPP6AXBqFGjNEu01K5dWwMu4P0aN26cvnfIHNOkSaPHccER+HwR6d+/vwZzvN/W+4Dn/eijj6Rly5bex+E1Wo+xJEmSRN9fq1T88ccfy82bN2X27Nl6H+CCq27dujJixAhJly5duD9HROTLY+C47SPNeEuUKHHfscKFC3u/xocvgpEVdK1jKD8DskAEIATbV155RT/ML126FK3nRsBDwLWyWwQ3KxgjA8PYpHUBAJMnT9bzRcBKmjSplm1PnTrlFwzxs3/88Yc3O6xTp45mwvDTTz9pRoqftW7I+JBNHz9+PNJzRdbqC5mk77kBvj98+LBm2uG9l7jYQIC13rtgIUPFBRAuLHxfAy6UcDyy8wX8jXyDJ14D/n5W0LVeA94PKysP7+cicuvWLbl69arfDceIyPzA6wniJm4PvL4fuhbfRhwEi8DGHBzDhzMg00QJEmPA+fPn1xI0srmoAhmgjIxsDcEK/8V4Lo4heCIQI3gg6MOXX34p7777rgYdlLrRCYwMGc0+llKlSmmmjMf+888/smjRIm+ZGTB22b59e/1Z64ZgjOfHzwX7PkVHZO9dsKyxV1zc+L4GjJd///33UZ7vg76G6P4cKgroAve94RgRmS3M5WO8sbKABoIJMiXcMP6IkjOCXrdu3SL9OWRSGF9Exla0aFHN3hB4UeZE1uzbbIQx6HLlynnLthCY5QECLTJdjAejWQsZr6V48eIa4HPlyvXQrxnNTTgnX/geJWdcjESHlUX6ZsiRQaUBDU7Hjh3zu6B4mNeACgAyaSu44jXgfXuQJioMSQT+zdF4NmG5/0UBEZnFY2AWG6fn8e7YsUMbfH744Qct+6JBCJ3K+FCP7jgvAqUVZFGaRXkS47MoRVswHovnQMf1oUOHdGwV5ehACEg//vijNh9h7Bgf/BaMsWLcGM1DyBSR6X7zzTf3NRNFB8Z7cY7odsb5oLkJ46PIyqMrbdq02k28YsUKOXfunFy5ciXKn8E4N7JIjE3jedHY9cknn8jYsWODfg14r9D1jLFhZM0Ye37rrbfktdde847vBgPvNbq+fW++7z8RmcnDUrO98OGK6TdoIkK216dPH23+qVWrVrR+HsEVGZ8VeJFtIRhbWbQFJWI0C7366qtSpkwZ+euvv/yyXwuy2dKlS2vXcmBWiKCOEjYCFqYUYRoPMnRkkcFC9ozGLpS1MecXv2fQoEF+jVVRSZAggQZQNEbhHOrVqxflz7Rp00YbnRBsUTHA+4es1ZqaFAyU8XEhg6lcKNPjQgXj9biAICKKrjCXl5pDPCZeTpCjjVr8/6eImahH/Upydvd2MV36EmXl6LwZYrqcjUNd8zpjSvPx/5vWGV1zujYWk3CTBCIispXH5flenAq8aJaKCLqgUe4lf3zPiMhpwhh4447IFv/PlCmTrefiFHzPiMhpPOJucSrwxsS0Hbfhe0ZE5CxxKvASEZH57j3gokCm4H68RERENmLGS0REtvK4fJCXgZeIiGzlcXnkZeAlIiJbhbk88HKMl4iIyEYMvERE5MpNEi5evKhr9GMPAezDjm1kre1UI7N9+3bdDx67tOFnsV8AtpaNLgZeIiJy5SYJzZs3l19//VX3iF+6dKlu4NOuXbsog27NmjWlevXqsnPnTt31DjvWYcOe6OIYLxERua65av/+/brFKgJnyZIl9djEiRN157zRo0dHuAvdO++8I126dJFevXp5jwW7HzkzXiIiitOl5lu3bsnVq1f9bjj2MJC5orxsBV2oVq2aZq7YNz4858+f1/uwN3q5cuV0H3Jstbply5agnpuBl4iI4nSpediwYZIsWTK/G449jLNnz2oADdzzPGXKlHpfeI4dO6b/HTBggLRt21YzZuy1jn3JDx8+HO3nZuAlIiJbeYL8X+/eveXKlSt+NxwLD0rAISEhkd4OHDjwQOcd9u9Sl+3bt5fWrVtLsWLFZNy4cVpqnjlzZrR/D8d4iYgoTkuUKJHeoqN79+7SqlWrSB+TI0cOSZ8+vZaOfd29e1c7nXFfeDJkyKD/zZ8/v9/xfPnyyalTpyS6GHiJiMiY5qo0adLoLSply5aVy5cvy+7du6VEiRJ6bN26dZrVlilTJtyfyZYtmzZdHTx40O/4oUOHpFatWtE+R5aaiYjIddOJ8uXLp9OCMFaLaUFbt27VaUFNmjTxdjSfPn1a8ubNq/cDytQ9evSQCRMmyPz58+XIkSPSt29fLV1jDnB0MeMlIiLXTSeCOXPmaLBFcxS6mRs2bKhB1XLnzh3Nbm/cuOE91rVrV7l586ZOK0JZukiRIjoPOGfOnBJdDLxERORKKVOmlLlz50Z4P0rL4V0koIHLdx5vsBh4iYjIVmFxJOONLQy8RERkK4+7466EeOJKsZ3IZlj5BpPwMR8wulMVnIiv0yxueZ0mY+Al18Kyc1gBB5PxscOIqfg6zeKW12kyTiciIiKyEQMvERGRjRh4iYiIbMTAS66FxpT+/fsb36DC12kWt7xOk7G5ioiIyEbMeImIiGzEwEtERGQjBl4iIiIbMfASERHZiIGXiIjIRgy8RERENmLgJVfDhtZERHZi4CXXCQsLk8GDB0umTJkkadKkcuzYMT3et29fmTFjhpjiwIEDEd63cuVKMcE///wjN27c8H5/8uRJGT9+vKxatUpMcefOHcmZM6fs378/tk+FYggDL7nOkCFDZNasWTJy5Eh57LHHvMcLFiwo06dPF1MUL15cJk+efN+Wcp07d5Z69eqJCfA6Zs+erV9fvnxZypQpI2PGjNHjU6ZMERMkTJiQlRnDMPCS6+CDetq0adK8eXOJHz++93iRIkUizRKdBhcX/fr1k9q1a8u5c+dk7969UqxYMVmzZo1s3rxZTPDjjz9KxYoV9ev58+dLunTpNOvF33jChAliijfffFNGjBghd+/eje1ToRiQICZ+CZGTnD59WnLlyhVuCRplPVM0btxYypUrJ61bt5YCBQrI9evXpVWrVpoRPvHEE2IClJmffPJJ/Rrl5QYNGki8ePHk2Wef1QBsil27dsnatWv1NRYqVEiSJEnid//ChQtj7dwoeAy85Dr58+fXjO/pp5/2O46MCRmhaW7fvi337t3TW4YMGSRx4sRiClxALV68WF5++WUdt37nnXf0+Pnz543aJD558uTSsGHD2D4NiiEMvOQ6KL+2bNlSM19kucgWDh48qOXJpUuXiim+/PJL6dixo5ZiDx06pKVmZL8IUJ999pnkyJFDTPhbNmvWTANu1apVpWzZsnocmaFJF1GffPJJbJ8CxSDuTkSuhIx30KBB8tNPP8m1a9e0EQkf4tWrVxdToBw5evRoDb6WS5cuSfv27WXFihVy9epVMcHZs2flzJkzOkaPMjPs3LlTM968efOKSS5cuKAXiZAnTx5JkyZNbJ8SPQAGXnIVNKe8//778sYbb0jmzJnFZPiAxodzeJDxvvbaa7afEz0YjM+/9dZbWpVBlQbQGPj666/LxIkTjRmzdwsGXnIdzN3dt2+fZMuWTdxwobFhwwY5evSolmTRiPTHH39oNoj3wemqVKkiISEhEd6/bt06MQGqFOhGnzRpkpQvX16PbdmyRbp06SIvvPCCMVOn3IJjvOQ6GAvcuHGj8YEXXb01a9aUU6dO6fxdfEAj8GJaCr6fOnWqOF3RokX9vkdXOsaycWGFcXxTLFiwQJv/Kleu7D2GaWKPP/64dq8z8DoLAy+5Tq1ataRXr17yyy+/SIkSJe6bmvHSSy+JCd5++20pWbKkjmOnSpXKexwdwG3bthUTjBs3LtzjAwYM0LF7U2DaFOYoB0qbNq3fyl3kDCw1k+tYDTjhQdkS025MgGC7bds2HedFposAjE7mEydO6JQqkz+wjxw5IqVLl5aLFy+KKVUa/D0xxmtNB8Nymcjq8RpRhibnYMZLrmM1p7jhdYZ3EfHf//7Xu+iEqbZv327UfGWsP41hAzQEonsbcCGF12jKuttuwsBLZChMjcIHNpbHtLJ5lF/79++v44MmwEpVvlDAw9SiH374QTe9MAVWqzp8+LDMmTPHu6xp06ZNddlTjPOSs7DUTK6D+buRwXxeEyCzrVGjhgYjfGhjvBf/TZ06tWzatEnHB50OC4IEDiNgbuvzzz9vzJxsNIxhPjIWd8mXL19snw7FAAZecp3AFY3wwXb8+HFJkCCBbr+GhfdNmk6EFax+/vln70IhzJKcB1tYYhyXgdcMDLxEIrqKEzYQQMcvF5aguAaLvmDZT2xbiQtEcjYGXqJ/YXpR3bp1tevXqZYsWRLtxzp12lTKlCk1CKFkniJFikgX0DClqxkXhNidCIuecHci5+OlE9G/rly5ojcnq1+/vt/3CEqB19ZWoHLqtCnM3bW6stE85gbcncgszHjJdQI3SLc6YbF+caVKlWTu3LliAowJvvfee1qmtHbtwTSbPn366DGsZEXOGKfHv0k0i6VPnz62T4diAAMvuU727Nkj7ITt3bu3MXNcCxYsqMtCVqhQ4b6dmdq1ayf79+8XJwpmVyVT9uTFJgj4ewXuIU3OxFIzuQ46mN0AGyOgRBkoWbJkjh7HxmuKbFzXl1PL6YGwCteePXsYeA3BwEuugy0BP/jgg/syW2vrtZkzZ4oJSpUqJd26ddMSurXO77lz56RHjx76Qe5U69ev936NCwisu42OdN9y+qeffirDhg0TU3Tq1Em6d++uc7PDW1+8cOHCsXZuFDyWmsl1sI8pxnQDF5D4888/dQwNY2qmrFeMblh0AGfJkkWP/f7775I7d25ZvHix5MqVS0xYw7hNmza6ipMvjIlixS5siWjq+uJW45xJ64u7BTNecg2MDeKDCre///7bby1ffHB99913RqzmZEFgxcIZq1ev9i4ziAUYqlWrFu1SbVyH7Da87Q2xShcCsincMjziFsx4yTWQNUQWcHDfwIED5T//+Y+t50UPDjsv1atXT0aOHOl3vGfPnvLNN9/IwYMHY+3ciCLCwEuusXHjRs120b2MjcWxEIPlscce08aVjBkzikmw6AJu58+fv29XJhPGslGlwPxWZPdlypTRYzt37tQ1qfE3NmUzCMBYPbJ7ZL/I9PHvFfOY0aWPiw9yDpaayTUwRxfwwYUxz8j25TUBsndsCIGya4YMGYwpL/tCYEWQnTJlind6FFYf69Chg3dc2wR4fdi8o2vXrjJ06FDvmC46vBF8GXidhRkvuRY2gj916pTcvn3byA5RBFuUYLn2tPPlz59fFz3BymToxsdevDly5JB9+/ZJ5cqVtTGQnIMZL7nOhQsXdDu55cuXh3u/KR2iuKAoV66cuIHpF1Go0gTuqgWJEiXSaXDkLGbX2ojCgXLd5cuXZceOHbo93ooVK3TeJ6bZBLPJQFyHrl5Tlr+M7CLqxRdf1CywQIECGpx8b6bAOO7evXvvO45/u9wq0HmY8ZLrrFu3TjteMfaJcV40qWDdYiwviEUX6tSpIya4efOmzmXFms3I/BImTOh3/9ixY8WkiyiUXBctWqSLhAwZMkTGjBkjpsBCKG+++ab+TTE6iAayL774Qv+9YqtAchYGXnIdlOas+brYVg5Z0zPPPKPbrf34449iCszhLVq0qH6NsUBfpjRaueUiCtULVGewwQXK6s2aNdMOfKzA1qRJk9g+PQoSAy+5cu4n5ndmy5ZNihQpIh999JF+jakaaEgyhe/SiqZyy0UUNG/eXG8IvNeuXQt3sZetW7fqRQjGfinu4hgvuc7bb7+tS0ZC//79tckqa9asul0gOkfJeRdRYF1EnT592riLqMCdiiJaYa1WrVr6+ilu43Qicj1kEFhSEcE3derU4mQNGjSQWbNmaakVX0dm4cKF4nSff/65rq2NTRJ2794tNWvWlIsXL+qCKHgfXn31VXET36lGFHex1Eyus2XLFr89apFBFC9eXEyALf+s8Vt8bboWLVp4v8auPSdPnjTmIorMxYyXXAfZUKZMmXRHG3xwY3ECN3PquOCdO3ckb968snTpUk6p+RczXmfgGC+5zh9//KF7m2Lt5oIFC2rn76hRo3SvUzdy6rggpkdheg2R0zDwkuugBNm5c2fN9I4ePSqvvPKKLqCBzmZsoOA2Ti56YW7riBEjjNlD+WGZMk3MdBzjJVfDikC9evXSjti+fftqFkzOsWvXLt19adWqVTqFKEmSJMY1kLnlIspNmPGSayHj7dSpk047wYIEKDsvW7Ystk+LgoDdebAtYI0aNXRBCTSU+d5MgUoMVugKdPXqVb8qzd9//83xXQdgcxW5DjLcr776Ssd6scoRFiXAtmrobnYjNzTkOLWBzIJVuc6ePXvf/F3ss4xGQTSakXOw1Eyus3nzZunRo4c0btyYU05cMi6IBjJsMuC0iwss+2n57bffNPj67qKFTRIQeMlZGHjJVZAZYLUjfBAz6P6PG4peTn2N6LjHhRFu4TX+Yf3miRMnxsq50YNjqZlcB2N/yH7QWGU6dPtu2LBBu7cxjo2yMkrsWNkqadKk4hZOLadjQRB8ROO8sSNRmjRp/Oajo/QcP378WD1HCh4zXnKd+vXry+LFi+Wdd94Rk+FDG0soYoP4W7du6Xg2AhCm3+B7rGdMcRt2W4KwsLDYPhWKQQy85DrY8H7QoEHacINlBgOnoHTp0kVM2QwCDUXI9FKlSuU9/vLLL0vbtm1j9dwoeIcPH9Ydp9BQFRiI+/XrF2vnRcFjqZlcJ7ISM8bSjh07JiZAsN22bZuOafuWWk+cOKHLZGJzCLdAad2JzVWWjz/+WDp27Kh9CenTp/driMPXpm2BaDpmvOQ6x48fFzdAVoTO10BYGhOB2E2cnl8MGTJEhg4dKu+9915snwrFAC6gQWSo6tWry/jx4/0yI2ygjj2Ia9euLSZwy8ISly5d0qVNyQwsNZPrvPHGG5HeP3PmTDEBMlus6IT/i2N8EOO9+C/KlZs2bYpwM3UnccvCEqGhoVKqVCnp0KFDbJ8KxQCWmsl1kD34wofzvn37NHMyaZOEzJkz67jul19+qQsxINvFBzhW6sL8Tydz28ISuXLl0rXEv//+e12TGjszmdgQ6BbMeIn+HQ9F80rOnDmlZ8+eYgJsmZc4cWIxETJdq8EovI8wa2GJqKobTuGWhkC3YOAl+tfBgwelcuXKcubMGTGlkxdTh1q0aCFVq1bVYGUKLixBTsZSM9G/sLqTSfu6Yo/huXPn6gYQWK3r1Vdf1SCMsV6nc/PCElau5IY1tk3FjJdcp1u3bn7f4/8CyHKxJWDLli1l0qRJYhJ09M6fP1+++OILWbdunWaJCMCmLLqACyZ0b+/fv1+/xxxlLB6CYQOTzJ49W0aNGqUNcvDMM8/oZh+vvfZabJ8aBYmBl1ynSpUqft+jBItSJRqrMCaYIIG5hSA0IqG5Cs1J4c3xdZqVK1fKSy+9pJsJlC9fXo9hRTI0lX377be6TKYJxo4dq81VnTt39r7OLVu2yOTJk3WOr+nLn5qGgZfIcGiyWrJkiZad0e2bLl06adq0qQwfPlycrlixYjplKvC1YM/lVatWGbOiE5qrBg4cKK+//vp9wwkDBgxwzaIwpmDgJdfBhxTGcrFmsy+U8DBNI1u2bGICZIMIttgQAll8o0aNNNt97rnnxBTo2v7ll1/u+1seOnRIChcurBcdprxOTHnDtKLAf7OYXmTK63QLc9ociaKpVatWuoZxoB07duh9pkBH8z///KNjg5jn+tFHHxkVdAFDBFiDORCOmbBAiAUBd968efcd/+qrr+676KC4z9zBLKII7NmzxztO5uvZZ5/VMTRTnDt3zvg1mbHLUrt27XQea7ly5bxjvNj6MLCJzslQZkZXOlYc8x3LXrt2bbgBmeI2Bl5yHUzDQKdvoCtXrji+4QhrFGP+LmAUCd9HxHqck6HhCBcXY8aMkd69e+uxjBkz6rinSas5NWzYUOcro8kKQweQL18+PYZxbnIWjvGS69StW1dXNsL0GmuRBQRcZBTXr1+X5cuXi1Ph9WBqFMqsvqs7+cL/5XHc6RcZgayLKdOyfCxp2r59e73IiGwFK3IOBl5yHUypwVhn8uTJpWLFinps8+bNmh1inmvBggXFqTZu3KilSDRT4evIVKpUSUyBTRGw8hjkzZvXbyUrE2ABFIxbM/CagYGXXOmPP/7QhTIw3xPZLzpgMb6bMmVKMcWpU6ckS5Ys92W9+L/877//LlmzZhUTstxOnTpp9cJaxQpZP6oXmOOKgGUCLOyCucqcr2sGBl6iCOADfdCgQbqNntPLzr7++usvPWZCqRkBFs1y2BChbNmyemz79u26chUCFXZmMgEWycA4NtbcLlGihCRJksTvfpPGs92AgZcokuYjlPecuoE6xnjR2RxYdsUGA1hWEePZTocAhPnKFSpU8DuOoYOaNWsa8RqBuxOZhV3NRBFw6jWpNY0GH8hoyHniiSe89yHLxXxlZIMmSJUqVbjlZBxLkSKFmNKhzpWpzMLAS2QYlF6tCwes6oRt8iz4ukiRIvLuu++KCfr06aMXGp999pmkT59ej2GxEGwegIsOJ8OFgzVUgHXEFy5cqA2B5HwsNRNFANNS0Hzl1FJz69at5YMPPjBivm5EMIf1yJEjcuvWLW+zGJrKEiVKdN+KTk5btxlZ+/fff6/zdSMaNiBnYsZLZChslRfe/sIXL17U6UYmBOT69euLqapVq6Y7aSHwWkuA+lYvfGEaHDkHAy+RoZo0aaKLhaA72xeWGMRuRd999504Xf/+/cVUn3/+ue4+hP2GMSe7QIECfuP15FwsNZMrNGjQQGbNmqVZHjYNwDQUlCMj07FjRxk8eLBjpxNhTjLW87UyJsuBAwd0kQ1MK3I6zEdGE1nmzJn1eyyhiB2Z0LWNNZxNgcx30aJFHOM1BHcnIldYunSpd2oJxj6xLnNUpkyZ4tigCxj3DK/UjCUIsWuRCZo1aybr16/3NlWhPIvg+5///EfnYJsCrxFB9/bt27pCV3h/V3IOlprJFbCMIBbRR+aAIg/KrRGNcQZuNu5UpUuXlmnTpuniEr6mTp2qizCYAHvU4nUC/qbYmxZZ/qpVq6RDhw7Sr18/MQEulLCyGkrP1n7DaPp76623JFOmTNKrV6/YPkUKAgMvuQKCDaadLFu2TEuTmIYS3gYCOGZK4MVqR8gA0ZmNFY8A28jt2rVLA5MJkL1bQwZr1qyRl156yXuhhak4pkBgxd9xw4YNujCIBX9f7MTEwOssHOMl18HUDJQlTdooPSJYeWvUqFH6X2tNamT+pmyeXqZMGa1i1KlTR6pXr67TbzBPGf9t1KiR/Pe//xUTPP3007rpPfaM9p3mhqlUxYsXj3T7R4p7mPGS62AVILfMh8QKVXPmzBFTYcN7TLPBxQU2EkDQBXRtWyVoE1y4cCHcC0X0LYRXuaG4jYGXXAfZw+XLl2XGjBmyf/9+PYYu2NDQUGN2swl08+ZNbczxZcI83sqVK8uff/6pGZ/vEpHoaPadeoNx35IlS0bZyR5X4dwxTIIxXbCC7fTp072bQ5BzsNRMrvPDDz9IjRo1tPRqZUUY90QDC8Y+UbozwY0bN6Rnz57adBTe1CETdidyy4YXW7ZskVq1akmLFi10Wlz79u11X+lt27bpHF9TmuXcgtOJyHWwpymacE6cOKHr3+KG8vOLL74oXbt2FVNgvWKsaIRpUcj0kB0NHDhQMmbMqHOZ3cTp+QV2X8KFA6YRoXMbF4goPWMLRAZd52HGS66DTBcbCaDz1RcyCJT0kCmaAGsXI8CiHIuMD2sV58qVSzcUwMbxJqxc5ZZ1t8kszHjJdRCEsJB+eKsg4QPaFFiT2Qo0eM343sqeNm3aFMtnR8HAtCGUmNm9bAYGXnIdLBeJRipMz0Cwxe3LL7+UNm3aSNOmTcUUCLrWPq7I7jHWC99++y2XHnQYrNOMaWDY+vCVV16Rb775RucwkzOx1Eyug+5ejH9iUQ1r6b2ECRPq2szDhw93bOdroHHjxkn8+PGlS5cuurgENkzA/93xgT127Fh5++23xS2c3lwFYWFh+nfEWtRYtxl/W8xVbt68uVSqVCm2T4+CwMBLroWxXOz8Ajlz5rxv5xcsvoBGJCy4YYKTJ0/K7t27dZwXC2m4iWljvJgehsrF0KFD5ZdffnFVh7oJGHiJDMySkNViaUFk9aasUhUeTAHDR5h10YSLC2SDmJeNlaxMhFXXMDSCbQPRMIcpcVipi5zDjEt5okfAydekKJ3//PPPYrp69ep5p0ZhURQsITlmzBg9jmlUpkBT1SeffCIvvPCCZMmSRV8bpsQdPnyYQdeBGHiJDIXFFrA6l8mQ8VWsWFG/nj9/vqRLl06zXgTjCRMmiCnwurDVYcGCBXXuLrYGxM5LGCIh5+GSkUSGQuPYzJkztSEHiywkSZLE7340WJkwTm9NAcOiEg0aNNAxeWwmgABsCqw9jR2mTOk3cDsGXiJDYa9aa/lL7N/qy5SF9dEotnjxYt0oYeXKlboqGZw/f96ItagtKDGTORh4iSLgxOCEcV2UI5EZrV+/XkyHcmuzZs004CIjtDYMQPZbrFgxMQlK6ZiLjcVfAje8QMmdnIN1CyKDmqsQbLBbD6AbO7zNEUyCeawIRNj4YsWKFd7jCMKYx2wKjFe3bt1ax3qx3Ck6mVOlSiXHjh3TzRPIWTidiFwPHaPYTCBPnjySL18+73GsaIV5vFiowCnwYYw1mNHdi6z33Llzrtl72GRYeax///66sprvnGRk/FgKdNKkSbF9ihQEBl5yncaNG8tzzz0nnTt31nmg2DwdOxXh/wqYH9mwYUNxKuxDi47eDBkyaCaYOXPmCC8ckC2ZsJDExIkTtayOcV2s7mRiCRbzlLF3NPaSxq5Eq1ev1n+3mE6ERjLTKxum4RgvuQ42CMDUDMBiCwi4mAP66aefypAhQxwdeKdNm6advUeOHNGlItu2bWvUxg+BsOY2xnNRckb51Ynj8tGBNZqR2SLwYtcpzN1F4MVa3MydnIeBl1znypUrkjJlSv0a44IItMgo6tSpo2s4Ox1WrAIsD4n1mE0OvEuXLtXSevny5cVkzz//vE4pwhg+xnrRTIZmK4xt40KLnIWBl1wHK/9gEQIEXwRelJfh0qVLkjhxYjEFVjoCZL9YkxrldexFjAzJlMwwU6ZMRl9Y+FYyrDL6m2++qWP527Zt09Wr2rdvH9unR0HiGC+5zocffqiZYNKkSbV0h3FANCJhrHDhwoXGTMNBaRJbyOH1INBiPBANOW+88YakSJFCl1Z0uuXLl2vHL9akxt+SyAk4nYhcp1OnTprxYlWnLVu2eFcDQlDCGK8punbtqms2o8nKd+cl7EfsO/XGyUqWLKkNVvjbIfNFFcP3ZgpUK9DBvHbtWn295GzMeIkMhYYcrOaEJhzfKSjoZsa2gNeuXROnq1atml5YoMkKc1wDS+gtW7YUE+CCEE2BKC9jKVBccFSuXFn34cX4duCWlhS3cYyXXAel1sggEzbB9evXw/1ARgk6UaJEYgIEIlQvcHFhsj59+uh/EXR37dolGzdulA0bNsjIkSO1YsMs2FlYaibXQROV7w3zP7GABsZ3Ma3IFNi1x9oyD5ANokEHH9ZVqlQRUxaWwFxst0C1Ahvfo3qB5UFRyeDKVc7DUjORiAakjh076jZrPXv2FBP8+uuvOg0FGyXgwgIdsDiGjHfr1q1GbCmHObwDBw6UoUOHSqFChXRM25cpGyVgPWpkubdu3dLxXpSYUWrGkIEpHepuwsBL9C/scYoPszNnzojT3blzR+fzDhs2TFc5QoaEMV0EYUxHwcpWJrAa4wKDjzVl6t69e2LK60ydOrUOk+BiqkKFChzXdTCO8RL9C3NdMYZmAmR+KEVi2pC1SpdpcHEBmEqEdbZNhiUhN2/erOO6vXv31uUjixYtqheKuFWvXj22T5GCwIyXXKdbt25+3+P/Ashyly1bpl2wpiw4j9WN0EQ1fPhwMRU2gECDVe7cucVNsCgKOp3nzJmjwySmZPZuwYyXXAfbqgWW8fABjgUloup4dhJk7+jQXrNmjZQoUUKSJEnid//YsWPF6Vq0aCEzZsww+uLCynitTmbcfvvtN0mePLnUrVtXx3vJWZjxEhkqss5ljH+i4crp3nrrLe3cRsZr6sUFYIcpjPGiU91qrEIzGTkTAy8ROZYbLi4A3egFChSI8nHoVsfiGqbM0zYVAy+5Arp5sdwemo2ww0tkUzBM2cOV3AfTp/bu3asrlFHcxTFecoV69ep5s4D69evH9ukQPRLMo5yBGS8RkSF81+SmuItLRhIREdmIpWZyBYztRndpPSypSET0qDDwkiuMHz/eb04kFh+oUaOGlC1bVo9hhxtsode3b99YPEuih8N1m52BY7zkOg0bNtRpKJ07d/Y7jhWrsNjE4sWLY+3ciB4Gx3idgYGXXCdp0qQ65SJXrlz3LcOH9W9N2CCezIKtD/FRbW2McPLkSVm0aJHkz5+f6zQ7EJuryHVSpUol33zzzX3HcQz3EcXF6XDW3srYM7pMmTK6xCmOT5kyJbZPj4LEMV5yHezf2qZNG13zFh9gsGPHDlmxYoV8/PHHsX16ROEu6jJu3Dj9ev78+ZIuXTpdc3zBggXSr18/3UuanIOBl1ynVatWki9fPpkwYYIsXLhQj+H7LVu2eAMxUVxy48YNHb+FVatWSYMGDXRzj2effVbLzuQsHOMlIorjChcurFWal19+WQoWLKjVGXTk7969W+rUqSNnz56N7VOkIHCMl1zt5s2bcvXqVb8bUVyDcvK7774r2bJl06qMNQ0O2S/WHidnYcZLrizb9ezZU+bNm6dzegNxU3GKi5DVnjlzRooUKaJlZti5c6dujJA3b97YPj0KAjNecp0ePXrodnHoBsXGCdOnT9eGq4wZM3o7R4nimvTp02t2awVdKF26NIOuAzHjJdfJmjWrBlhsJo5sAR2jmNP72WefyRdffCHfffddbJ8i0X1DIhMnTpT169fL+fPnJSwszO9+bmXpLOxqJtfBWszWyj4IvNbazBUqVOC0DIqTQkNDdTy3UaNGmuVyaUhnY+Al10HQPX78uGa+KNNhrBcfZt9++60kT548tk+P6D5Lly7VSkz58uVj+1QoBnCMl1yndevWup4t9OrVSyZPniyJEyeWd955R8d/ieKaTJkyeefxkvNxjJdc5c6dO1KzZk2ZOnWq5M6dW49hAQLMh8Q4L+ZLEsU1y5cv1wVf8O/26aefju3ToYfEUjO5SsKECeXnn3/2O4YPMn6YUVxWsmRJbbDCMAk2SsC/Y1/cQ9pZmPGS66CkjGlEw4cPj+1TIYqWatWqyalTp7TJCus0BzZXtWzZMtbOjYLHjJdc5+7duzJz5kzde7dEiRKSJEkSv/vHjh0ba+dGFJ5t27bJ9u3bdfEMcj4GXnKdffv2SfHixfXrQ4cO+d3HaRoUF6H7HnvykhlYaiYiiuMwhxerqw0dOlQKFSp03xgv5qOTczDwEhHFcdYykYEVGXx84xjXF3cWlpqJiOL4FDjAVKI8efLE9ulQDGDGS0QUx6VJk0YbrKy55+RsXLmKiCiOa9GihcyYMSO2T4NiCEvNRERxHKfAmYWBl4gojuMUOLNwjJeIiMhGHOMlIiKyEQMvERGRjRh4iYiIbMTAS0REZCMGXiIiIhsx8BIREdmIgZeIiEjs8/8AOWIXz/u7FQkAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 500x300 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# This cell is based on https://seaborn.pydata.org/examples/many_pairwise_correlations.html\n",
    "\n",
    "# Compute the correlation matrix\n",
    "corr = df.corr()\n",
    "\n",
    "# Set up the matplotlib figure\n",
    "f, ax = plt.subplots(figsize=(5, 3))\n",
    "\n",
    "# Generate a custom diverging colormap\n",
    "cmap = sns.diverging_palette(230, 20, as_cmap=True)\n",
    "\n",
    "# Draw the heatmap with the mask and correct aspect ratio\n",
    "sns.heatmap(corr, cmap=cmap, vmax=0.3, center=0, square=True, linewidths=0.5)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several relationships we can identify:\n",
    "\n",
    "- The RMS spot size positively correlates to the RMS wavefront error, as expected.\n",
    "- There is no correlation between refractive index and radius of curvature, as expected.\n",
    "- There is a strong negatively correlation between the RMS spot radius and the radius of curvature.\n",
    "\n",
    "\n",
    "Note that these relationships are not universally true, but are based on the Monte Carlo input distributions used, the system tested, etc. With different inputs, or a different optical system, the results would undoubtedly differ."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
